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Houston,  Texas 


Based  on  the  support  provided  by  the  aforementioned  grant, 
considerable  advances  were  made  in  several  structural  mechanics 
problems  with  stochasticity. 

A  new  method  for  the  solution  of  problems  involving  material 
variability  was  developed.  The  material  variability  was  modeled 
as  a  stochastic  process.  The  Karhunen-Loeve  expansion  of  random 
processes  was  used  to  represent  the  material  variability  in  a 
computationally  expeditious  manner.  The  well-known  deterministic 
finite  element  method  has  been  employed  to  discretize  the 
differential  equation  which  governs  the  nodal  response  random 
variables.  An  related  spectral  expansion  of  these  random 
variables  was  adopted  in  terms  of  the  basis  in  the  space  of 
second  nodal  random  variables.  This  method  yielded  a 
representation  of  the  response  surface  in  terms  of  the  polynomial 
chaos.  The  coefficient  in  this  representation  was  such  that  it 
involved  enough  information  about  the  process  so  that  one  could 
reproduce  its  probability  distribution  function.  The  method  has 
been  applied  to  a  plain-stress  problem  which  involves  a  curved 
geometrical  boundary.  The  representation  of  the  random  field 
over  the  curved  domain  was  accomplished  by  solving  the  related 
integral  equation  using  a  Galerkin  formulation.  Interestingly, 
the  result  of  the  representation  is  independent  of  the  mass  size 
which  was  employed,  and  converged  quite  rapidly  as  the  number  of 
terms  in  the  Karhunen-Loeve  expansion  increased.  Even  more 
encouraging  was  the  fact  that  the  analytical  results  were  found 
in  extremely  good  agreement  with  data  produced  by  a  Monte-Carlo 
study  of  the  problem.  The  findings  of  this  research  effort  have 
been  summarized  in  the  paper: 

"Stochastic  Finite  Element  Analysis  With  Curved 
Boundaries,"  by  R.  G.  Ghanem  and  P.  D.  Spanos, 

Proceedings  of  the  Sixth  International  Conference  on 
Application  of  Statistics  and  probability  in  Civil 
Engineering,  Mexico  City,  Mexico,  1991,  June  7-9,  1991, 
pp.  158-165. 

An  effort  was  undertaken  to  develop  a  formulation  of 
stochastic  concepts  towards  a  unification  of  various  finite 
element  techniques.  Specifically,  methods  for  the  solution  of 
differential  equations  with  random  processes  and  coefficients 
have  been  addressed.  The  method  which  was  advocated  treats  the 
random  aspects  of  the  problem  as  an  added  dimension.  In  this 


manner,  a  formulation  for  the  stochastic  finite  element  method 
has  been  derived  which  can  be  construed  as  a  natural  extension  of 
the  deterministic  finite  element  method.  Finite  element 
representation  al^png  the  random  dimension  was  achieved  by  two 
spectral  expansions.  One  of  them  was  used  to  represent  the 
coefficients  of  the  differential  equation  which  modeled  the 
random  material  property,  and  the  other  one  was  used  to  represent 
the  random  solution  process.  The  new  concepts  were  implemented 
by  using  a  number  of  computational  models  as  simple  engineering 
systems.  The  new  method  indeed  involves  generalizing  the 
concepts  of  finite  element  approximation  to  abstract  spaces  of 
which  the  usual  Euclidian  space  is  a  special  case.  The 
deterministic  case  can  be  regarded  as  a  digression  of  this 
formalism  to  the  particular  case  where  the  space  of  elementary 
events  involves  only  a  single  element,  and  where  the  probability 
density  function  induced  on  the  associated  a-algebra  is  the 
uniform  distribution.  The  findings  of  these  research  efforts 
have  been  summarized  in  the  paper: 

"A  Spectral  Formulation  of  Stochastic  Finite  Elements," 

by  R.  G.  C-hanem  and  P.  D.  Spanos,  Proceedings  of  the 

Tenth  International  Invitational  Unification  of  Finite 

Element  Methods  Symposium,  July  18-19,  1991,  Worcester 

Polytechnic  Institute,  Worcester,  MA,  pp.  59-82. 

Another  numerical  method  for  dealing  with  problems  of 
stochastic  mechanics  has  been  pursued  and  it  involves  a  small 
number  of  random  parameters.  This  method  is  analogous  to  the 
Monte-Carlo  simulation  method,  but  it  is  more  efficient.  In 
fact,  the  method  treats  a  stochastic  problem  as  an  ensemble  of 
deterministic  ones.  After  solving  a  number  of  deterministic 
problems,  statistical  analysis  is  performed  to  deduce  the 
necessary  parameters  characterizing  the  random  nature  of  the 
solution.  Of  course,  this  method  is  often  the  only  option  which 
is  available  to  solve  complicated  stochastic  problems.  However, 
indiscriminate  use  of  the  method  is  not  advocated  due  to  the 
significant  computational  cost.  The  new  method  has  been  used  for 
the  determination  of  the  eigenvalues  of  a  beam  bending  problem 
with  random  parameters.  The  beam  is  assumed  to  be  clamped- 
clamped  of  unit  length,  and  its  rigidity  is  a  truncated  normal 
process  with  mean  equal  to  one  and  with  exponential  auto¬ 
correlation  function.  This  kind  of  problem  is  quite  difficult  to 
be  treated  either  analytically  or  by  other  numerical  methods.  In 
fact,  the  available  algorithms  can  be  applied  in  the  case  of  very 
small  randomness  and  they  are  quite  costly  computationally.  The 
new  method  has  been  used  to  determine  the  first  three  eigenvalues 
of  this  problem.  It  was  found  that  the  analytical  results  are  in 
very  good  agreement  with  the  results  obtained  by  numerical 
simulations  of  this  problem.  Further,  it  has  been  found  that  the 
new  method  can  be  applied  to  a  wide  class  of  problems  dealing 
with  random  variables  and  stochastic  processes.  The  findings  of 
this  research  effort  have  been  summarized  in  the  paper: 


"Pseudo-Simulation  Method  For  Stochastic  Problems,"  by 
B.  A.  Zeldin  and  P.  D.  Spanos,  Proceedings  of  the  Sixth 
Probabilistic  Mechanics,  Structural  and  Geotechnical 
Reliability  Specialty  Conference,  American  Society  of 
Civil  Engineers,  Boulder,  CO,  July  7-9,  1992,  pp.  37- 
40  . 

Another  approach  to  dealing  with  problems  of  the  numerical 
solution  of  stochastic  mechanics  problems  has  been  pursued. 
Specifically,  the  Finite  Difference  Method  for  discretization  of 
stochastic  continuous  media  has  been  addressed.  The  Neumann 
expansion  and  perturbation  methods  used  for  solving  the  systems 
of  the  associated  algebraic  equations  have  been  studied. 

Further,  work  has  been  done  into  the  dependence  on  the  mass  size 
of  the  discretization.  It  has  been  found  that  a  mixed 
formulation  which  involves  both  strain  and  stresses  as 
independent  variables  improves  the  performance  of  the  method  and 
reduces  the  dependence  on  the  mass  size.  The  method  has  been 
used  to  study  the  behavior  of  a  beam  which  is  subjected  to 
deterministic  load  and  involves  bending  rigidity  which  is  a 
normal  random  process.  It  has  been  found  that  the  mixed 
formulation  includes  the  convergence  of  the  Neumann  expansion, 
and  it  minimizes  its  dependence  on  the  mass  size.  This  is 
particularly  true  with  regards  to  the  variability  of  the 
displacement,  even  when  the  coefficient  of  variation  of  the 
bending  rigidity  is  quite  large.  The  findings  of  this  research 
effort  have  been  summarized  in  the  paper: 

"Stochastic  Mixed  Finite  Difference  Method,"  by  P.  D. 

Spanos  and  B.  A.  Zeldin,  Proceedings  of  the  Sixth 
Probabilistic  Mechanics,  Structural  and  Geotechnical 
Reliability  Specialty  Conference,  American  Society  of 
Civil  Engineers,  Boulder,  CO,  July  7-9,  1992,  pp.  804- 
807. 

A  class  of  Galerkin-type  projection  procedures  applied  for 
random  media  for  stochastic  mechanics  problems  has  been 
developed.  The  use  of  a  space  of  simple  random  variables  as 
projective  has  been  adopted.  Further,  it  has  been  shown  that  in 
this  case  this  procedure  is  equivalent  to  a  generalization  of 
sampling  technique.  Then,  the  optimality  of  this  indirect 
sampling  has  been  investigated.  It  has  been  shown  that  this 
method  is  a  generalized  stratified  sampling  method.  That  is,  the 
approximation  of  the  problem  from  the  space  of  sample  random 
variables  produces  an  error  with  mean  and  conditional 
expectation,  given  sigma-algebra  induced  by  this  partition,  equal 
zero.  It  has  been  shown  that  this  error  has  zero  estimate  from 
the  set  of  indicator  functions  of  given  stratification.  Some 
stochastic  mechanics  problems  have  been  considered  utilizing  the 
proposed  method.  It  has  been  shown  that  the  Loeve-Karhunen 
expansion  is  a  versatile  tool  for  the  approximation  of  a 
stochastic  field  through  a  set  of  random  variables.  This 


expansion  is  optimal  and,  therefore,  a  small  number  of  random 
variables  can  approximate  them  adequately.  Several  examples  have 
demonstrated  that  the  proposed  method  can  be  applied  for  treating 
a  wide  class  of  problems  dealing  with  random  variables  and 
stochastic  processes.  The  findings  of  this  research  effort  have 
been  summarized  in  the  paper: 

"Indirect  Sampling  Method  For  Stochastic  Mechanics 
Problems,"  by  P.  D.  Spanos  and  B.  A.  Zeldin, 

Proceedings  of  the  Sixth  International  Conference  of 
Structural  Safety  and  Reliability,  Innsbruck,  Austria, 

August  8-13,  1993  (to  appear). 

In  terms  of  personnel,  it  is  stated  that  in  addition  to  the 
supervising  professor,  two  research  assistants  were  involved  in 
the  projects.  One  was  supported  exclusively  by  the  research 
grant  from  AFOSR,  and  another  was  jointly  supported  by  this  AFOSR 
grant  and  discretionary  funds  from  Rice  University.  The  two 
research  assistants  and  the  supervising  professor  formed  a  very 
coherent  group  which  worked  with  enthusiasm  in  this  challenging 
technical  field. 
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Based  on  the  support  provided  by  the  aforementioned  grant, 
considerable  advances  were  made  in  several  structural  mechanics 
problems  with  stochasticity. 

A  new  method  for  the  solution  of  problems  involving  material 
variability  was  developed.  The  material  variability  was  modeled 
as  a  stochastic  process.  The  Karhunen-Loeve  expansion  of  random 
processes  was  used  to  represent  the  material  variability  in  a 
computationally  expeditious  manner.  The  well-known  deterministic 
finite  element  method  has  been  employed  to  discretize  the 
differential  equation  which  governs  the  nodal  response  random 
variables.  An  related  spectral  expansion  of  these  random 
variables  was  adopted  in  terns  of  the  basis  in  the  space  of 
second  nodal  random  variables.  This  method  yielded  a 
representation  of  the  response  surface  in  terms  of  the  polynomial 
chaos.  The  coefficient  in  this  representation  was  such  that  it 
involved  enough  information  about  the  process  so  that  one  could 
reproduce  its  probability  distribution  function.  The  method  has 
been  applied  to  a  plain-stress  problem  which  involves  a  curved 
geometrical  boundary.  The  representation  of  the  random  field 
over  the  curved  domain  was  accomplished  by  solving  the  related 
integral  equation  using  a  Galerkin  formulation.  Interestingly, 
the  result  of  the  representation  is  independent  of  the  mass  size 
which  was  employed,  and  converged  quite  rapidly  as  the  number  of 
terms  in  the  Karhunen-Loeve  expansion  increased.  Even  more 
encouraging  was  the  fact  that  the  analytical  results  were  found 
in  extremely  good  agreement  with  data  produced  by  a  Monte-Carlo 
study  of  the  problem.  The  findings  of  this  research  effort  have 
been  summarized  in  the  paper: 

1}  "Stochastic  Finite  Element  Analysis  With  Curved 
Boundaries,"  by  R.  G.  Ghanem  and  P.  D.  Spanos, 

Proceedings  of  the  Sixth  International  Conference  on 
.  Application  of  Statistics  and  probability  in  Civil 
Engineering,  Mexico  City,  Mexico,  1991,  June  7-9,  1991, 
pp.  158-165. 

An  effort  was  undertaken  to  develop  a  formulation  of 
stochastic  concepts  towards  a  unification  of  various  finite 
element  techniques.  Specifically,  methods  for  the  solution  of 
differential  equations  with  random  processes  and  coefficients 
have  been  addressed.  The  method  which  was  advocated  treats  the 
random  aspects  of  the  problem  as  an  added  dimension.  In  this 


manner,  a  formulation  for  the  stochastic  finite  element  method 
has  been  derived  which  can  be  construed  as  a  natural  extension  of 
the  deterministic  finite  element  method.  Finite  element 
representation  along  the  random  dimension  was  achieved  by  twc 
spectral  expansions.  One  of  them  was  used  to  represent  the 
coefficients  of  the  differential  equation  which  modeled  the 
random  material  property,  and  the  other  one  was  used  to  represent 
the  random  solution  process.  The  new  concepts  were  implemented 
by  using  a  number  of  computational  models  as  simple  engineering 
systems.  The  new  method  indeed  involves  generalizing  the 
concepts  of  finite  element  approximation  to  abstract  spaces  of 
which  the  usual  Euclidian  space  is  a  special  case.  The 
deterministic  case  can  be  regarded  as  a  digression  of  this 
formalism  to  the  particular  case  where  the  space  of  elementary 
events  involves  only  a  single  element,  and  where  the  probability 
density  function  induced  on  the  associated  a-algebra  is  the 
uniform  distribution.  The  findings  of  these  research  efforts 
have  been  summarized  in  the  paper: 

2)  "A  Spectral  Formulation  of  Stochastic  Finite  Elements," 
by  R.  G.  Ghanem  and  P.  D.  Spanos,  Proceedings  of  the 
Tenth  International  Invitational  Unification  of  Finite 
Element  Methods  Symposium,  July  18-19,  1991,  Worcester 
Polytechnic  Institute,  Worcester,  MA,  pp.  59-82. 

Another  numerical  method  for  dealing  with  problems  of 
stochastic  mechanics  has  been  pursued  and  it  involves  a  small 
number  of  random  parameters.  This  method  is  analogous  to  the 
Monte-Carlo  simulation  method,  but  it  is  more  efficient.  In 
fact,  the  method  treats  a  stochastic  problem  as  an  ensemble  of 
deterministic  ones.  After  solving  a  number  of  deterministic 
problems,  statistical  analysis  is  performed  to  deduce  the 
necessary  parameters  characterizing  the  random  nature  of  the 
solution.  Of  course,  this  method  is  often  the  only  option  which 
is  available  to  solve  complicated  stochastic  problems.  However, 
indiscriminate  use  of  the  method  is  not  advocated  due  to  the 
significant  computational  cost.  The  new  method  has  been  used  for 
the  determination  of  the  eigenvalues  of  a  beam  bending  problem 
with  random  parameters.  The  beam  is  assumed  to  be  clamped- 
clamped  of  unit  length,  and  its  rigidity  is  a  truncated  normal 
process  with  mean  equal  to  one  and  with  exponential  auto¬ 
correlation  function.  This  kind  of  problem  is  quite  difficult  to 
be  treated  either  analytically  or  by  other  numerical  methods.  In 
fact,  the  available  algorithms  can  be  applied  in  the  case  of  very 
small  randomness  and  they  are  quite  costly  computationally.  The 
new  method  has  been  used  to  determine  the  first  three  eigenvalues 
of  this  problem.  It  was  found  that  the  analytical  results  are  in 
very  good  agreement  with  the  results  obtained  by  numerical 
simulations  of  this  problem.  Further,  it  has  been  found  that  the 
new  method  can  be  applied  to  a  wide  class  of  problems  dealing 
with  random  variables  and  stochastic  processes.  The  findings  of 
this  research  effort  have  been  summarized  in  the  paper: 
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3) 


"Pseudo-Simulation  Method  For  Stochastic  Problems,"  by 
B.  A.  Zeldin  and  P.  D.  Spanos,  Proceedings  of  the  Sixth 
Probabilistic  Mechanics,  Structural  and  Geotechnical 
Reliability  Specialty  Conference,  American  Society  of 
Civil  Engineers,  Boulder,  CO,  July  7-9,  1932,  pp.  37- 
40. 

Another  approach  to  dealing  with  problems  of  the  numerical 
solution  of  stochastic  mechanics  problems  has  been  pursued. 
Specifically,  the  Finite  Difference  Method  for  discretization  of 
stochastic  continuous  media  has  been  addressed.  The  Neumann 
expansion  and  perturbation  methods  used  for  solving  the  systems 
of  the  associated  algebraic  equations  have  been  studied. 

Further,  work  has  been  done  into  the  dependence  on  the  mass  size 
of  the  discretization.  It  has  been  found  that  a  mixed 
formulation  which  involves  both  strain  and  stresses  as 
independent  variables  improves  the  performance  of  the  method  and 
reduces  the  dependence  on  the  mass  size.  The  method  has  been 
used  to  study  the  behavior  of  a  beam  which  is  subjected  to 
deterministic  load  and  involves  bending  rigidity  which  is  a 
normal  random  process.  It  has  been  found  that  the  mixed 
formulation  includes  the  convergence  of  the  Neumann  expansion, 
and  it  minimizes  its  dependence  on  the  mass  size.  This  is 
particularly  true  with  regards  to  the  variability  of  the 
displacement,  even  when  the  coefficient  of  variation  of  the 
bending  rigidity  is  quite  large.  The  findings  of  this  research 
effort  have  been  summarized  in  the  paper: 

4)  "Stochastic  Mixed  Finite  Difference  Method,"  by  P.  D. 

Spanos  and  B.  A.  Zeldin,  Proceedings  of  the  Sixth 
Probabilistic  Mechanics,  Structural  and  Geotechnical 
Reliability  Specialty  Conference,  American  Society  of 
Civil  Engineers,  Boulder,  CO,  July  7-9,  1992,  pp.  804- 
807. 

A  class  of  Galerkin-type  projection  procedures  applied  for 
random  media  for  stochastic  mechanics  problems  has  been 
developed.  The  use  of  a  space  of  simple  random  variables  as 
projective  has  been  adopted.  Further,  it  has  been  shown  that  in 
this  case  this  procedure  is  equivalent  to  a  generalization  of 
sampling  technique.  Then,  the  optimality  of  this  indirect 
sampling  Has  been  investigated.  It  has  been  shown  that  this 
method  is  a  generalized  stratified  sampling  method.  That  is,  the 
approximation  of  the  problem  from  the  space  of  sample  random 
variables  produces  an  error  with  mean  and  conditional 
expectation,  given  sigma-algebra  induced  by  this  partition,  equal 
zero.  It  has  been  shown  that  this  error  has  zero  estimate  from 
the  set  of  indicator  functions  of  given  stratification.  Some 
stochastic  mechanics  problems  have  been  considered  utilizing  the 
proposed  method.  It  has  been  shown  that  the  Loeve-Karhunen 
expansion  is  a  versatile  tool  for  the  approximation  of  a 
stochastic  field  through  a  set  of  random  variables.  This 
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expansion  is  optimal  and,  therefore,  a  small  number  of  random 
variables  can  approximate  them  adequately.  Several  examples  have 
demonstrated  that  the  proposed  method  can  be  applied  for  treating 
a  wide  class  of  problems  dealing  with  random  variables  and 
stochastic  processes.  The  findings  of  this  research  effort  have 
been  summarized  in  the  paper: 

5)  "Indirect  Sampling  Method  For  Stochastic  Mechanics 
Problems,"  by  P.  D.  Spanos  and  B.  A.  Zeldin, 

Proceedings  of  the  Sixth  International  Conference  of 
Structural  Safety  and  Reliability,  Innsbruck,  Austria, 

August  8-13,  1993  (to  appear). 
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1)  Stochastic  Finite  Element.  Analysis  with  Curved  Boundaries 

R.G.  Ghancin  *  P  D.  Spanos  1 


Abstract 

An  original  method  for  the  solution  of  problems  involving  material  variability  it  proposed.  The 
material  properly  is  modeled  as  a  slorh.u»Uc  proe  u.  The  K arhunen-Loeve  expansion  is  used  to 
represent  this  profcM  in  a  computationally  expedient  manner.  The  standard  deterministic  finite 
clniiriH  motli>.u|  is  i-mploynl  to  >li«rr«lt/.4<  the  flilT.-rout  mi  l tons  governing  ths  nodal  rrjponse 

random  variables.  A  spectral  expansion  of  these  nimlom  variables  ia  then  adopted  in  terms  of 
a  busts  ill  the  space  of  second  order  random  variables.  The  method  yields  an  expression  for  the 
response  surface  in  terms  of  the  Polynomial  Chaoses.  The  coefficients  in  the  expansion  ara  such 
that  they  involve  enough  information  about  the  response  process  so  as  to  be  possible  to  reproduce 
its  probability  distribution  function.  The  method  ia  applied  to  a  plane-streas  problem  involving 
a  rnrvo’d  g«>r«m>uni-nl  boundary .  lb  pr-M.  Hl  .it  »«  u  of  ihr  raurlom  field  over  the  curved  domain  is 
acLOmphshed  by  solving  the  related  integral  equation  using  u  O'aicfkm  formuiaiion.  1  ha  rnmuU'utg 
representation  ia  independent  of  the  inesh  sir.-;  employed  and  converges  rapidly  as  the  number  of 
terms  in  the  Kurhunrit- Locve  expansion  i6  r. .creased. 


1  Introduction 

The  analysis  of  engineering  systems  with  uncertain  properties  has  witnessed  a  considerable  resurgence 
in  recent  years,  in  particular  with  relation  to  systems  involving  many  degrees  of  freedom  and  requiring 
recourse  to  the  finite  element  method  for  their  analysis  and  design.  A  comprehensive  treatment 
of  this  class  of  problems  can  be  accomplished  by  breaking  down  the  complexity  into  two  separate 
issues.  The  first  one  consists  of  adequately  representing  the  uncertainty  in  the  system  properties 
for  implementation  within  a  computational  framework.  The  Karhuncn-Locve  expansion  is  used  in 
this  paper  as  an  optimal  such  representation.  It  effectively  replaces  the  random  process  by  a  set 
of  uncorrelated  random  variable*  while  delegating  the  corresponding  spatial  dependence  to  a  set  of 
deterministic  functions.  The  second  issue  involved  in  the  solution  process  is  obtaining  a  representation 
for  the  response  process.  This  task  is  accomplished  by  first  identifying  a  complete  basis  for  such  a 
representation,  and  then  by  formulating  the  problem  in  such  a  way  that  the  coefficients  of  this  basis 
in  the  representation  can  be  numerically  computed.  In  this  paper,  earlier  formulations  of  the  ideas 
presented  above  (Spanos  and  Gbanem,  1989,  1990,  1991)  arc  extended  to  deal  with  situations  where  the 
geometry  of  the  domain  under  consideration  is  irregular,  involv.ng  curved  boundaries.  The  efficiency 
of  the  proposed  method  in  treating  such  problems  is  demonstrated  by  its  application  to  the  analysis  of 
a  curved  thin  plate.  In  addressing  this  problem,  it  is  reminded  that  the  ultimate  goal  of  a  stochastic 
finite  clement  analysis  is  the  Calculation  of  certain  statistics  of  the  response  procoai.  Tltcsa  statistics 
can  be  in  the  form  of  either  statistical  moments,  or  probability  distribution  function,  or  some  other 
measure  of  the  reliability  of  the  system.  A*  a  first  step  in  the  solution  procedure,  the  variational 
formulation  of  the  finite  element  method  is  used  to  obtain  a  spatially  discrete  form  of  the  problem. 
Following  that,  the  Polynomial  Chaos  expansion  is  used  to  derive  a  representation  of  the  response 
process.  Statistical  moment*  and  probability  distribution  functions  can  then  obtained. 
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Figure  1:  Plato  with  Random  Rigidity,  Exponential  Covariance  Model. 

2  Finite  Element  Formulation 

Consider  th  ■  thin  plate  shown  in  Figure  (1  j.  Its  modulus  of  elasticity  is  assumed  to  be  the  realization  of 
a  two-dimensional  Gaussian  random  process  with  known  mean  value  £  and  known  covariance  function 
C(Xi,Xj).  Further,  it  is  assumed  that  the  external  excitation  is  deterministic  and  of  unit  magnitude. 
Let  the  domain  A  of  the  plate  be  discretized  into  N  four-noded  quadrilateral  fir.i-e  elements,  each 
element  having  eight  degrees  of  freedom.  The  strain  energy  y  stored  in  each  element  A *  can  be 
expressed  as 

v  =  i  /  <rT(x)  «(x)  <IA'  .  (1) 

2  J,\< 

where  dA‘  is  a  diferential  element  in  /l*.  Further,  <r(x)  and  e(x)  denote  the  stress  and  the  strain 
vectors  respectively,  ,  s  a  function  of  the  location  x  within  each  element.  Assuming  linear  elastic 
material  behavior,  the  stress  may  be  expressed  in  terms  of  t!.e  strain  as 

<r  =  D'  *  (2) 

where  D*  is  the  matrix  of  constitutive  relations.  Here  tr  and  «  -.re  the  vectors  of  stress  and  strain  a* 
given  by  the  equation 

VT  -  (  "r,  Cr,  ’•r.r,  I  (3) 

*T  =  [  fx.  (rj  Tex,  1  .  (4) 

where  a,%  is  the  stress  along  direction  x,  and  rXi  is  the  strain  along  that  same  direction.  For  the  plane 
stress  problem  considered  herein,  D°  is  given  by  the  equation 

pktyr\  1  Pe  0 

D'  =  l‘<  !  0  =  F.'(x)  P'  ,  (5) 

1  - "«  0  n 

where  P*  it  a  deterministic  matrix,  is  the  elemental  Poisson  ratio,  and  £'(x)  is  the  elemental 
modulus  of  elasticity.  The  two  dimensional  displacement  vector  u(x)  representing  the  longitudinal  and 
transverse  displacements  within  each  element  may  bo  expressed  in  terms  of  the  nodal  displacements 
of  the  element  in  the  form 

u(sc)  *  H*(r,,r,)  U*  ,  (6) 

where rj)  is  the  local  interpolation  matrix,  U*  is  the  random  nodal  response  vector, and  rI  and 
rj  are  local  coordinates  over  the  element.  Substituting  equations  (5)  and  (0)  into  equation  (1)  gives 

V  -  \  /  £'(x)  eT(x)  P'  efx)  d/t*  . 

L  J  A* 


(7) 
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The  strain  within  an  element  is  related  to  the  displacements,  longitudinal  end  transverse,  through  the 
relation 


«(x) 


t)i, 


0 


0 


d 

dzt 


u(x)  . 


11 

l  dz j  dx\  J 


(8) 


Using  equation  (6),  equation  (7)  is  rewritten  as 

e(x)  =  D'  U*  , 


(0) 


where  matrix  B  involve*  derivatives  of  the  interpolation  matrix  II.  SubiUiuiing  equation  (0)  back 
into  equation  (7)  and  performing  a  coordinate  transformation,  leads  to 


1"  =  \vT  ['  ('  £(r,.r,)  D'T(--.,r3)  P'  D*(r,,r,)  jJ'Idr,  dr}  V  , 

l  Jo  Jo 


(10) 


where  |J'|  denotes  the  determinant  of  the  Jacobian  of  the  transformation  that  maps  an  arbitrary 
clement  (e)  onto  the  four  noded  square  with  sides  equal  to  one.  The  total  strain  energy  V  is  obtained 
by  summing  the  contributions  from^all  the  elements.  This  procedure  gives 


V  =  l  £  U‘T  /'  /’  £(r..'-l)B*T(rl,r,)P*B*(r„r,)|J‘|rfr1  dr,  U* 
2  ft  Jo  Jo 
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The  local  representation  of  the  response  is  related  to  the  global  representation  through  the  following 
transformation 

U*  =  C*  U  .  (II) 

where  C*  is  a  rectangular  permutation  matrix  of  zeros  and  ones  reflecting  the  connectivity  of  the 
elements  and  the  topology  of  the  mesh.  Using  equation  (12),  the  following  expression  for  the  total 
energy  stored  in  the  system  is  obtained 

V  =  i  UT  K  U  .  (13) 

.  2 

Before  evaluating  the  integrals  in  equation  (11),  the  random  process  representing  the  modulus  of 
elasticity  of  the  plate  must  be  replaced  by  a  representation  which  somehow  decouples  its  spatial 
dependence  from  its  random  dependence.  The  Karhunen-Loeve  expansion  (Spanos  and  Ghanem, 
1988)  is  used  herein  as  an  optimal  such  representation.  Accordingly,  £(r, ,r3)  is  replaced  by 


M 

£(ri.,‘3,  =  53  \AT?t/*(ri.ri)  .  (H) 

i  »  1 

where  A*  and  arc,  respectively,  the  eigenvalues  and  eigenfunctions  associated  with  the 

covariance  function  of  the  random  process.  The  next  stage  in  the  computations,  therefore,  involves 
solving  l lie  integral  eigenvalue  problem  associated  with  the  covariance  kernel.  That  is, 

/n(x,,y,)  =  LC{  ri,yi;*j.yj)  fJn.v j)  dxa  dj/3  •  05) 

The  kernel  used  in  this  paper  is  defined  by  the  equation 

=  *-  I*'-"1/*'  *  >*>-«<'*»  ,  (16) 


where  and  fr?  arc  the  correlation  distances  in  the  i\  and  Xj  directions  respectively. 
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3  Numerical  Solution  of  the  Integral  Equation 

Subdividing  the  domain  A  of  the  plate  into  N  finite  clement*  A* ,  equation  (IS)  become* 

*  , 

K  t.Vt)  =  ]C  /  /»(*i.Va)  dA*  •  (17) 

«-i  -',4* 

Interpolating  for  the  value  of  the  unknown  function  within  an  clement  in  term*  of  it*  nodal  value* 
results  in  the  following  expression 


/.(!.»)  =  ll'trt.r,)  i;  ,  (18) 

where  H‘(ri,rj)  is  the  interpolation  matrix  in  terms  of  local  coordinate*  r|  and  r  j,  and  f*  is  the  vector 
of  nodal  value*  for  the  unknown  function  associated  with  clement  (e).  For  this  particular  problem, 
bilinear  interpolation  is  used  over  four-noded  quadrilateral  clement.  The  matrix  II*(ri,rj)  It  then 
given  by  the  equation 

H*(rltr,)  =  i  1  (1  -  r,)(l  -  r,)  (1  +  r,)(l  -  r,)  (19) 

(l+r,)(l  +  r,)  (l-r,)(l  +  r,))  . 

Substituting  equation  (18)  and  performing  a  transformation  from  global  to  local  coordinates,  equation 
(17)  becomes 

s  r 

K  /„(*!.».)  =  5Z  /  C(*,.V,;*J.Vj)H'(r„r,)|J*|  <fA*  t*  ,  (20) 

sal 

where  [J*j  is  the  Jacobian  of  the  coordinate  transformation.  A  system  of  algebraic  equation*  i*  obtained 
from  equation  (28)  by  requiring  the  corresponding  error  to  be  orthogonal  to  all  the  interpolation 
functions  used.  That  is, 

CD  =  A  B  D  ,  (21) 

where  now  the  j,k  column  of  D  is  the  j"1  eigenfunction  calculated  at  the  nodal  points  and 

A,;  =  f,j  Xi  .  (22) 

Matrices  C  and  B  are  obtained  by  assembling  matrices  Ce/  and  B ,/  where 

c'{  =  f  jA/  C(i„!q;x1,yJ)  II*T(r,.r,)  H'(r„r,)  rfA*  dA>  ,  (23) 

and 

n,/  =  J  C(jr,.y,;r,,yi)  !rr(r,.r,)  H/(r,,r,)  dA*  ,  (24) 

with  (zi,yi)  and  (rt,ri)  denoting  the  global  and  local  coordinates  of  a  point  in  A*  respectively.  The 
assembly  procedure  just  mentioned  consists  of  combining  entries  corresponding  to  the  same  node 
(Akin,  1982). 

4  Spectral  Stochastic  Finite  Elements 

The  Karhuncn-Locvc  expansion  for  the  modulus  of  elasticity  may  be  substituted  into  equation  (11) 
to  transform  equation  (13)  into 

1  « 

v  =  5 ur  12  f* K(l’ u 

2  t»i 


(25) 
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The  integration!  involved  in  equation  (25)  may  be  performed  either  analytically  or  using  some  numer¬ 
ical  quadrature  scheme.  The  work  performed  by  the  externally  applied  forces  is 


N  , 

V‘  =  Y  u'  r(x)  f(x)  ,IA'  . 

A'  r 

-  UT  y  CtT  /  ir  T(x)  f(x)  dA‘  =  Urf. 

S  =  !  J 

(2G) 

where 

N 

f  =  Y'  CrT  /  ll'(x)  f(x)  dA‘  . 

,=i  J  '• 

(27) 

Minimizing  the  total  potential  energy  ( V  -  t")  with  respect  to  U,  leads  to  the  equation 


K'°>  +  £  f*  K'l> 


U 


(28) 


The  re'ponse  vector  U  is  expanded  along  the  basis  given  by  the  Polynomial  Chaoses  (Ghanem  and 
Spanos,  1990)  as, 

U  =  £c,  *,({£,}].  (29) 

i=0 

Substituting  equation  (29)  into  equation  (28),  and  requiring  the  ettor  resulting  from  truncating  the 
scries  at  the  PiK  term  to  be  orthogonal  to  the  P  +  1  Polynomial  Chaoses  ♦.[{fr)]^o>  results  in  a 
system  of  linear  algebraic  equations  of  the  form 


p 

E 

t*0 


M 


£  <{**.[{M)*r({6}]>  K“> 


c,  =  <*,({{/})*>.  j  0...P  , 


(30) 


which  can  be  solved  for  the  coefficients  c, 


5  Numerical  Results 

A  curved  plate  is  shown  in  Figure  (1).  The  curved  side  is  a  ninety  degree  arc  of  a  circle  of  unit 
radius.  The  length  of  the  straight  edges  is  equal  to  1.25.  The  standard  deviation  of  the  longitudinal 
displacement  at  node  A,  using  two  terms  in  the  K-L  expansion  for  the  material  stochastic! ty,  is 
shown  in  Figure  (2).  It  is  plotted  against  the  standard  deviation  of  the  modulus  of  elasticity.  The 
results  corresponding  to  four  terms  in  the  K-L  expansion  are  shown  in  Figure  (3).  Note  the  excellent 
convergence.  The  probability  distribution  functions  corresponding  to  one  of  the  response  variable  at 
node  A  arc  depicted  in. Figures  (4)-(7).  Results  corresponding  to  two  and  four  terms  in  the  Karliuncn- 
Locve  expansion  and  up  to  third  order  Polynomial  Chaos  are  shown. 

The  two  dimensional  process  representing  the  modulus  of  elasticity  of  the  plate  is  simulated  In 
such  a  way  as  to  accommodate  the  non-uniform  spatial  distribution  of  the  nodal  points.  The  issue 
is  addressed  by  using  the  Karhunen-Locve  expansion  to  simulate  a  truly  continuous  random  Held. 
To  this  end,  the  eigenvalues  and  eigenfunctions  of  the  covariance  kernel  are  computed  as  described 
in  the  previous  section.  The  random  field  is  then  simulated  using  equation  (14)  with  the  number 
of  terms  equal  to  the  number  of  nodes  in  the  system.  The  orthogonal  random  variables  appearing 
in  that  equation  are  obtained  as  pseudorandom  computer  generated  uncorrelatod  variates,  with  xero 
mean  and  unit  variance.  The  resulting  simulated  random  field  is  not  as  sensitive  to  the  mesh  slxe  and 
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Figure  2:  Normalized  Standard  Deviation  of  Longitudinal  Displacement  at  Corner  /t  of  the  Curved  Plate,  vereue 
Standard  Deviation  of  the  Modulus  of  Elasticity;  Two  Terms  in  the  K-L  Expansion;  Exponential  Covariance; 

ffm.e  =  104- 


Figure  3i  Normalized  Standard  Deviation  ofLongitudinal  Displacement  at  Corner  A  of  the  Curved  Plate,  versus 
Standard  Deviation  of  the  Modulus  of  Elasticity;  Four  Terms  in  the  K-L  Expansion;  Exponential  Covariance; 
r«u  =  10.4. 
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Figure  d:  Longitudinal  Displacement  at  the  Free  End  of  the  Curved  Plate;  Probability  Density  Function 
Using  30,000-Samples  MSC,  and  Using  Third  Order  Homogeneous  Chaos;  Two  Terms  in  the  K-L  Expansion; 
Exponential  Covariance. 
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Figure  St  Longitudinal  Outplacement  at  the  Free  End  of  the  Curved  Plate;  Tail  of  the  Probability  Density 
Function  Using  30,000-Sample  MSC,  and  Using  Third  Order  Homogeneous  Chaos;  Two  Terms  in  the  K-L 
Expansion;  Exponential  Covariance. 
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Figure  Oi  Longitudinal  Displacement  at  the  Free  End  of  the  Curved  Plate;  Probability  Density  Function 
Using  30,000-Sample  MSC,  and  Using  Third  Order  Homogeneous  Chaos;  Four  Terms  in  the  K-L  Expansion; 
Exponentisl  Covarisnce. 
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Figure  7:  Longitudinal  Displacement  at  tile  Free  End  of  the  Curved  Plate;  Tail  of  the  Probability  Density 
Function  Using  30,000-Sample  MSC,  and  Using  Third  Order  Homogeneous  Chaos;  Four  Terms  in  tbs  K-L 
Expansion;  Exponential  Covariance. 
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nodal  point  distribution  as  the  field  obtained  using  more  conventional  procedures.  The  results  from 
using  the  Monte  Carlo  simulation  method  are  superimposed  on  the  same  plot  as  the  analytical  results. 
Observe  the  good  agreement  between  the  analytical  and  the  simulated  results  even  for  large  values  of 
the  coefficient  of  variation. 
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Abstract 

A  method  Tor  the  solution  of  differential  equations  with  random  processes 
as  coefficients  is  developed.  The  method  relies  on  viewing  the  random  aspect 
of  the  problem  as  an  added  dimension,  and  on  treating  random  variables  and 
processes  as  functions  defined  over  that  dimension.  This  way,  a  formulation  for 
the  stochastic  finite  element  method  is  developed  which  is  a  natural  extension 
of  the  deterministic  finite  element  method.  Finite  element  representation 
along  the  random  dimension  is  achieved  via  two  spectral  expansions.  One  of 
them  is  used  to  represent  the  coefficients  of  the  differential  equation  which 
model  the  random  material  properties,  the  other  is  used  to  represent  the 
random  solution  process. 
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1.  INTRODUCTION 

Although  the  first  quantitative  ideas  of  probability  theory  took  form  as  far  back 
as  the  seventeenth  century,  at  the  gambling  tables,  with  Pascal  and  Fermat,  it  was 
not  until  well  into  this  century  that  the  body  of  knowledge  known  as  probability 
became  a  mathematical  discipline.  In  the  interim,  some  of  the  greatest  scientific 
minds,  including  Gauss,  Laplace,  and  Poincare,  as  well  as  many  others,  were  laying 
the  foundation  for  the  theory  by  compiling  observations  and  related  theorems  on  the 
chance  occurrence  of  events.  Several  attempts  to  baptize  the  accumulated  ideas  as 
a  branch  of  mathematics  were  vehemently  opposed  by  prominent  mathematicians, 
as  lacking  the  rigorous  foundation  fit  of  a  mathematical  theory.  However,  the 
postulation  of  the  uncertainty  principle  early  in  this  century  created  an  urgency 
for  providing  a  sound  mathematical  framework  for  probablistic  concepts.  By  this 
time,  important  contributions  in  that  direction  had  been  made  by,  among  others, 
Poincare  and  Borel.  But  it  was  with  Kolmogorov,  in  his  Foundation  of  the  Theory 
of  Probability  (1933),  that  an  axiomatic  foundation  of  the  theory  was  presented 
and  that  the  subject  matter  finally  gained  universal  acceptance  as  a  branch  of 
mathematics. 

The  usefulness  of  this  axiomization  cannot  be  overestimated.  Indeed,  it  provided 
the  connection  between  probability  as  a  collection  of  observations  of  natural  phe¬ 
nomena  and  mathematical  reasoning,  thus  providing  a  whole  new  set  of  perspectives 
»nd  tools  with  which  to  view  and  approach  related  problems.  Specifically,  as 
related  to  the  development  of  stochastic  finite  elements,  the  most  significant  aspect 
of  mathematical  probability  is  the  association  of  random  variables,  which  are  the 
elementary  ingredients  of  the  theory,  with  functions  defined  over  topological  spaces. 
Once  this  association  has  been  established,  the  well  developed  field  of  functional 
•°aIysis  could  be  used  in  analyzing  and  operating  on  these  random  variables.  This 
connection  with  functional  analysis  already  carries  the  ingredients  for  a  unification 
°f  stochastic  finite  elements  with  deterministic  finite  elements.  Indeed,  a  major  part 
°f  the  modern  development  in  the  theory  of  finite  elements  draws  intimately  from 
functional  analysis,  so  that  at  a  certain  level  of  abstraction,  both  the  deterministic 
*^d  the  stochastic  finite  element  methods  have  the  same  theoretical  foundation. 

,  ce  this  unification  is  established,  the  deterministic  finite  element  method  can  be 
v,e'w«d  a,  a  jp^iaj  digression  of  the  stochastic  finite  element  method,  whereby 
*orne  0f  (he  functional  spaces  involved  have  a  particular  structure.  Unlike  the 
e  erministic  case,  however,  where  functions  are  usually  defined  with  respect  to 
1  tebesgue  measure,  when  dealing  with  random  entities  a  more  general  concept  of 
Measure  is  called  for.  Whereas  a  Lebesgue  measure  coincides,  usually,  with  the  more 
mtuitive  notion  of  differential  volume,  probability  measures  are  abstractions  of  such 
0  '”I,es-  The  practical  efTect  of  this  difference  is  that  whereas  in  the  deterministic 
Ijf3*  *  discretization  of  a  function  with  respect  to  its  natural  measure  induces  a 
.  *cretiaation  of  the  physical  space  with  respect  to  this  same  measure,  and  thus 

uc«  a  finite  element  mesh,  a  similar  discretization  in  the  probabilistic  case,  with 
re*pect  to  the  probability  measure,  does  not  carry  a  parallel  physical  consequence. 
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Although  at  first  this  may  be  seem  a  weak  point  in  the  development  of  a  finite; 
element  theory,  it  is  quickly  reminded  that  recent  efforts  in  deterministic  fiait*, 
elements  have  been  advocating  similar  abstract  discretizations,  the  p-method  ani 
various  spectral  methods  being  cases  in  point.  In  this  respect,  the  spectral  theory  o(, 
the  stochastic  finite  element  method,  as  developed  in  the  sequel,  can  be  viewed  as  t 
natural  extension  of  these  developments  to  the  case  of  random  operators.  In  the  next, 
section,  a  coherent  mathematical  framework  is  presented  which  is  a  natural  setting 
for  the  analysis  of  random  operator  equations.  Next,  the  theory  of  representation  of 
stochastic  processes  is  expanded  with  special  emphasis  on  two  spectral  expansions,' 
namely  the  Karhunen-Loeve  and  the  Polynomial  Chaos  expansions.  These  are  then, 
used  in  the  following  section  to  develop  the  stochastic  finite  element  method. 

2.  THE  MATHEMATICAL  MODEL 

The  class  of  problems  dealt  with  in  this  study  is  not  of  the  conventional  engineering 
kind  in  that  it  involves  concepts  of  a  rather  abstract  and  mathematical  nature.  It  is 
both  necessary  and  instructive  to  introduce  at  this  point  the  mathematical  concepts 
which  are  used  in  the  sequel. 

The  Hilbert  space  of  functions  (Oden,  1979)  defined  over  a  domain  D,  with  values 
on  the  real  line  R,  is  denoted  by  H.  Let  (O,  'I',  P)  denote  a  probability  space.  By 
that  is  meant  that  fl  is  a  space  of  elementary  events,  ♦  is  the  o- field  generated 
by  fl,  or  loosely  speaking,  the  space  consisting  of  the  various  combinations  of  the 
elements  of  fl,  and  finally,  P  is  the  probability  measure  defined  on  <£.  Let  x  be  an 
element  of  D  and  9  be  an  element  of  fl.  Then,  the  space  of  functions  mapping  fl 
onto  the  real  line  is  denoted  by  0.  Each  map  ft  —  R  defines  a  random  variable. 

The  inner  products  over  H  and  over  0  are  defined  using  the  Lebesgue  measure 
and  the  probability  measure,  respectively.  That  is,  for  any  two  elements  A,(x)  and 
fij(x)  in  H,  their  inner  product  (  h,(x)  ,  A,(x)  )  is  defined  as 

(Ai(x),A,(x))  =  Johi{x)h,(x)dx  .  (1) 

The  domain  D  represents  the  physical  space  over  which  the  problem  is  defined. 
Similarly,  given  any  two  elements  a(0)  and  0(9)  in  0,  their  inner  product  is  defined 
as 

(«{<?), /?(*))  =  Jna(9)0(O)dP  (2) 

where  dP  is  a  probability  measure.  Under  very  general  conditions,  the  integral 
in  equation  (2)  is  equivalent  to  the  average  of  the  integrand  with  respect  to  the 
probability  measure  dP,  so  that 

(a(0), /?(<>))  =  <a(0)m>  (3) 

where  <.>  denotes  the  operation  of  mathematical  expectation.  Any  two  elements 
of  the  Hilbert  spaces  defined  above  ?re  said  to  be  orthogonal  if  their  inner  product 
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^wishes.  A  random  process  may  then  be  described  as  a  function  defined  on  the 
product  space  D  x  fl.  Viewed  from  this  perspective  a  random  process  can  be  treated 
w  a  curve  in  either  of  H  or  ©. 

The  physical  model  under  consideration  involves  a  medium  whose  properties 
exhibit  random  spatial  fluctuations  and  which  is  subjected  to  a  random  external 
excitation.  The  mathematical  representation  of  this  problem  involves  an  operator 
equation 

A(x,0)(u(x,fl)]  =  /(x,0)  (4) 

where  A(x,  0)[.]  is  some  operator  defined  on  H  x 0.  In  other  words,  A  is  a  differential 
operator  with  coefficients  exhibiting  random  fluctuations  with  respect  to  one  or  more 
of  the  independent  variables.  The  aim  then  is  to  solve  for  the  response  u(x,B)  as 
a  function  of  both  its  arguments.  With  no  loss  of  generality,  A  i  assumed  to  be  a 
differential  operator  whose  random  coefficients  are  restricted  -eing  second  order 
random  processes.  This  is  not  a  severe  restriction  for  practical  problems,  since  most 
physically  measurable  processes  are  of  the  second  order  type.  Then,  each  one  of 
these  coefficients  at(x,  9)  can  be  decomposed  into  a  purely  deterministic  component 
and  a  purely  random  component  in  the  form 

ak{x,9)  =  a*(x)  +  atk(x,9)  (5) 

where  ak(x)  is  equal  to  the  mathematical  expectation  of  the  process  at(x,9),  and 
at(x,9)  is  a  zero-mean  random  process,  having  the  same  covariance  function  as  the 
process  at(x,9).  Equation  (4)  can  then  be  written  as 

(L(x)  +  n(x,*))[u(x,*)l  =  f(x,9),  (6) 

where  L(x)[.]  is  a  deterministic  differential  operator  and  II(x,0)[.)  is  a  differential 
operator  whose  coefficients  are  zero-mean  random  processes.  Before  a  solution  to 
equation  (6)  is  sought,  it  is  essential  to  clarify  what  is  meant  by  such  a  solution. 

“  will  prove  instructive  to  start  with  the  deterministic  finite  element  method 
^d  see  how  the  related  concepts  can  be  generalized.  A  finite  element  solution  to  a 
deterministic  problem  governed  by  a  certain  differential  equation  consists  basically 
of  computing  the  value  of  the  dependent  variables  on  a  discrete  mesh  induced  in  the 
•pace  spanned  by  the  independent  variables.  This  is  probably  the  most  widespread 
mterpretation  of  a  finite  element  solution;  it  has  been  crucial  in  disseminating  the 
Method  as  a  powerful  analysis  and  design  tool  into  engineering  practice.  An  alterna¬ 
te  viewpoint  which  will  prove  to  be  more  amenable  to  the  required  generalizations, 
u  that  a  solution  to  a  finite  element  pro  >lem  consists  in  evaluating  the  value  of  the 
efficients  in  the  expansion  of  the  solution  along  a  certain  basis  in  an  appropriate 
unetional  space.  The  finite  element  procedure  will  consist  in  choosing  a  suitable 
asis  and  then  computing  optimal  values  of  the  coefficients  with  respect  to  this  basis. 

perspective,  the  finite  element  mesh  is  naturally  induced  with  specific 
°>ces  of  these  bases.  With  other  choices,  however,  the  expansion  coefficients  do 
D°t  necessarily  carry  an  obvious  physical  interpretation.  In  the  stochastic  case, 
°D*  of  the  independent  variables  spans  the  space  of  elementary  events,  which  can 
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only  be  discretized  with  respect  to  a  probability  measure,  the  result  lacking  m 
intuitive  appeal.  In  this  case,  the  appeal  of  the  second  interpretation  of  a  fiajfi 
element  solution  is  obvious.  The  problem  then  becomes  one  of  identifying  a  suiUlS 
basis  in  the  .space  H  x  ©  over  which  the  solution  is  defined,  and  of  determina^ 
a  meaningful  optimality  criterion  for  computing  the  coefficients  in  the  associate 
expansion.  Obviously,  the  basis  functions  in  this  case  will  be  random.  By  simulating 
realizations  of  these  functions,  corresponding  realizations  of  the  solution  process  ag, 
be  obtained.  Alternatively,  by  defining  a  suitable  inner  product  over  the  space  ^ 
random  variables,  various  statistics  or,  equivalently,  norms  of  the  solution  procai 
may  be  evaluated. 

3.  REPRESENTATION  OF  STOCHASTIC  PROCESSES 

Similarly  to  the  case  of  the  deterministic  finite  element  method,  whereby  functions 
are  represented  by  a  denumerable  set  of  parameters  consisting  of  the  values  of  the 
function  and  its  derivatives  at  the  nodal  points,  the  problem  encountered  in  tie 
stochastic  case  is  that  of  representing  a  random  process  by  a  denumerable  set  J 
random  variables,  thereby  discretizing  the  process. 

In  the  deterministic  case  discretization  of  the  domain  has  a  physical  appeal 
The  domain  in  the  stochastic  case  does  not,  however,  have  a  physical  meaning  that 
permits  a  sensible  discretization.  In  this  context  the  functional  analysis  foundatioa 
of  the  finite  element  method  becomes  useful  as  it  can  be  extended  to  deal  with! 
random  functions.  Two  of  the  most  useful  expansions  for  random  processes  are  tit 
Karhunen-Loeve  expansion,  and  the  Polynomial  Chaos  expansion.  The  first  requires 
knowledge  of  the  covariance  structure  of  the  process  under  consideration,  while  tb< 
second  one  is  more  general.  The  difference  between  these  two  expansions  can  be 
loosely  compared  to  that  between  a  modal  expansion  and  a  Fourier-type  expansion 
of  a  system  response.  Although  the  former  has  better  convergence  properties,  the 
latter  is  more  general  and  does  not  require  knowledge  of  the  properties  of  the  system- 
These  two  expansions  are  discussed  next. 

3.1  Karhunen-Loeve  Expansion 

The  major  conceptual  difficulty  from  the  viewpoint  of  the  class  of  problems  con¬ 
sidered  herein,  involves  the  treatment  of  functions  defined  on  these  abstract  spaces, 
namely  random  variables  defined  on  the  or— field  of  random  events.  The  most  widely 
used  method,  the  Monte  Carlo  simulation,  consists  of  sampling  these  functions  *• 
randomly  chosen  elements  of  this  <r-field,  in  a  random,  collocation-like,  scheme- 
Obviously,  a  quite  large  number  of  points  must  be  sampled  if  a  good  approximat'00 
is  to  be  achieved.  Alternatively,  these  functions  could  be  expanded  in  a  Fourier-typ*' 
series  as 

w(x.0)  =  £  /-(x)  , 

n*sl 


(7) 
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•here  {(,(0)}  '*  4  set  rar,dom  variables  to  be  determined,  A„  is  some  constant, 
and  (/.{*)}  >*  an  orthonormal  set  of  deterministic  functions.  This  is  exactly  what 
the  K&rhunen-Loeve  expansion  achieves.  The  expansion  was  derived  independently 
by  *  number  of  investigators  (Karhunen,  1947;  Loeve,  1948;  Kac  and  S;egert,  1947). 

Let  w(x,9)  be  a  random  process,  function  of  the  position  vector  x  defined  over 
the  domain  D,  with  9  belonging  to  the  space  of  random  events  ft.  Let  u>(x) 
denote  the  expected  value  of  to(x,  9)  over  all  possible  realizations  of  the  process,  and 
C(xi,xj)  denote  its  covariance  function.  By  definition  of  the  covariance  function,  it 
it  bounded,  symmetric  and  positive  definite.  Thus,  it  has  the  spectral  decomposition 
(Courant  and  Hilbert,  1953) 

C(x„x,)  =  £  A.  /„(x.) /.(x,)  (8) 

ns  1 


where  A*  and  /„(x)  are  the  eigenvalue  and  the  normalized  eigenvector  of  the  covari¬ 
ance  kernel,  respectively.  That  is,  they  are  the  solution  to  the  integral  equation 

C(x,,X2)  /„(x)  dx,  =  A„  /„(xi)  .  (9) 

Due  to  the  symmetry  and  the  positive  definiteness  of  the  covariance  kernel  (Loeve, 
1977),  its  eigenfunctions  are  orthogonal  and  form  a  complete  set.  They  have  further 
been  normalized  so  that  the  following  equation  holds, 

[  /„(x)  /m<x)  dx  =  Snm  ,  (10) 

J  D 

where  S„„  is  the  Kronecker  delta.  Clearly,  w(x,")  can  be  written  as 

w{x,6)  =  u<(x)  +  o(x,0),  (11) 

•'he.e  o(x,0)  is  a  process  with  zero  mean  and  covariance  function  C(Xi,Xj).  The 
Process  o(x,0)  can  be  exp^.ided  in  terms  of  the  eigenfunctions  /„(x)  as 

O(x,0)  =  £  UB)  >/X,;«(x).  (12) 

1*1 

^*<*nd  order  properties  of  the  random  variables  can  be  determined  by  multiplying 
lh  sides  of  equation  (12)  by  o(xj  ,  0)  and  taking  the  expectation  on  both  sides. 
Pacifically,  it  is  found  that 

C(x„Xl)  =  <a(x,,0)  a(x,,0)> 

=  £  £  <U<>)  U(0)  >  TaTaT  /n(x.)  /m(x3)  . 

nsl  ms| 


(13) 
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Then,  multiplying  both  sides  of  equation  (14)  by  /*(xj),  integrating  over  the  dom*^ 
D,  and  making  use  of  the  orthogonality  of  the  eigenfunctions,  yields 

J  C(x,,x2)  /*(x3)  dxt  =  A*  /*(x,)  (U) 

=  E  <*.(*)  &(«)  >  j  *0  ■ 

n=l 

Multiplying  once  more  by  fi(xt)  and  integrating  over  D,  gives 

J  h  (x.)  ft  (x,)  dx,  =  £  E<U8)  M*)  >  i/UT  («J 

U  »=tl 

Then,  using  equation  (10)  leads  to 

At  &kt  —  \j^k  A|  <£k{0)  (i(0)  >  .  (16) 

Equation  (16)  can  be  rearranged  to  give 

<tk{9)  U0)  >  =  Su  .  (17) 

Thus,  the  random  process  w(x,9)  can  be  written  as 

u>(x,0)  =  uS(x)  +  £  (n(9)  \f\n  f„(x)  .  (18) 

nscl 

where, 

<U9)  >  =  0  ,  <U9)  U{9)>  =  *.»  ,  («) 

and  A„,  /„(x)  are  solution  to  equation  (9).  Truncating  the  series  in  equation  (18) 
at  the  Mtk  term,  gives 

w(x,6)  =  w(x)  +  E  U0)  VA„  /n(x)  .  (20J1 

n=0 

An  explicit  expression  for  {„(#)  can  oe  obtained  by  multiplying  equation  (12)  ty 
/„(x)  and  integrating  over  the  domain  D.  That  is, 

U*)  =  ~  lDo(x,9)  fn{x)  dx  .  {M 

It  is  well  known  from  functional  analysis  that  the  steeper  a  bilinear  form  dec*# 
to  zero  as  a  function  of  one  of  its  arguments,  the  more  terms  a  t  needed  in  '*4' 
spectral  representation  in  order  to  reach  a  preset  accuracy.  Noting  that  the  Fouri* 
transform  operator  is  a  spectral  representation,  it  may  be  concluded  that  the  f*»^ 


the  autocorrelation  function  tends  to  zero,  the  broader  is  the  corresponding  spectral 
density,  and  the  greater  the  number  of  requisite  terms  to  represent  the  underlying 
random  process  by  the  Karhunen-Loeve  expansion. 

For  the  special  case  of  a  random  process  possessing  a  rational  spectrum,  the 
integral  eigenvalue  problem  can  be  replaced  by  an  equivalent  differential  equation 
that  is  easier  to  solve  (Van  Trees,  1968).  In  the  same  context,  it  is  reminded  that  a 
necessary  and  sufficient  condition  for  a  process  to  have  a  finite  dimensional  Markov 
realization  is  that  its  spectrum  be  rational  (Kree  and  Soize,  1986).  Further,  note 
that  analytical  solutions  for  the  integral  equation  (10)  are  obtainable  for  some  quite 
important  and  practical  forms  of  the  kernel  C(xi,x3)  (Juncosa,  1945;  Slepian  and 
Poliak,  1961;  Van  Trees,  1968).  In  the  general  case,  however,  the  integral  equation 
must  be  solved  numerically.  Various  techniques  are  available  to  this  end  (Ghanem 
and  Spanos,  1991). 

3.2  Homogeneous  Chaos 

It  is  clear  from  the  preceding  discussion  that  the  implementation  of  the  h'arhunen- 
Loeve  expansion  requires  knowicdgeof  the  covariance  function  of  the  process  being 
expanded.  As  far  as  the  system  under  consideration  is  concerned,  this  implies  that 
the  expansion  can  be  used  for  the  random  coefficients  in  the  operator  equation. 
However,  it  cannot  be  implemented  for  the  solution  process,  since  its  covariance 
function  and  therefore  the  corresponding  eigenfunctions  are  not  known.  An  alterna¬ 
tive  expansion  is  clearly  needed  which  circumvents  this  problem.  Such  an  expansion 
«>uld  involve  a  basis  of  known  random  functions  with  deterministic  coefficients  to  be 
found  by  minimizing  some  norm  of  the  error  resulting  from  a  finite  representation. 
This  should  be  construed  as  similar  to  the  Fourier  series  solution  of  deterministic 
differential  equations,  whereby  the  series  coefficients  are  determined  so  as  to  satisfy 
*°me  optimality  criterion.  To  clarify  this  important  idea  further,  a  general  functional 
form  of  the  solution  process  is  written  as 

«*  =  h  [f.(#),*]  (22) 

*kere  h|.J  js  4  nonlinear  functional  of  its  arguments.  In  equation  (22),  the  random 
P'ocesses  involved  have  all  been  replaced  by  their  corresponding  Karhunen-Loeve 
^Presentations.  It  is  clear  now  that  what  is  required  is  a  nonlinear  expansion 
0  A!-l  in  terms  of  the  set  of  random  variables  i,(9).  If  the  processes  defining 
.  *  °Perator  are  Gaussian,  this  set  is  a  sampled  derivative  of  the  Wiener  process 
oob,  1953).  In  this  case,  equation  (22)  involves  functionals  of  the  Brownian 
Motion.  This  is  exactly  what  the  concept  of  Homogeneous  Chaos  provides.  This 
y'jceP1  was  first  introduced  by  Wiener  (1938)  and  consists  of  an  extension  of 
.  erra’s  work  on  the  generalization  of  Taylor  series  to  functionals  (  Volterra, 

2  ).  Wiener’s  contributions  were  the  result  of  his  investigations  of  nonlinear 
liQCl'°na  *  brownian  motion.  Based  on  Wiener’s  ideas,  Cameron  and  Martin 
' )  constructed  an  orthogonal  basis  for  nonlinear  functionals  in  terms  of  Fourier- 
crinite  functionals. 
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3.2.1  Definitions  and  Properties 

Let  {{1(0)}“,  be  a  set  of  orthonormal  Gaussian  random  variables.  Consider  the 
space  r,  of  all  polynomials  in  {6(0)}“i  of  degree  not  exceeding?.  Let  repreaeat 
the  set  of  all  polynomials  in  tr  orthogonal  to  r,_t.  Finally,  let  f,  be  the  spaa 
spanned  by  T?.  Then,  the  subspace  f,  of  ©  is  called  the  p,h  Homogeneous  Chao, 
and  F,  is  called  the  Polynomial  Chaos  of  order  p. 

Based  on  the  above  definitions,  the  Polynomial  Chaoses  of  any  order  p  consist 
of  all  orthogonal  polynomials  of  order  p  involving  any  combination  of  the  random 
variables  {4. (^)) “ x -  It  is  clear,  then,  that  the  number  of  Polynomial  Chaoses  of 
order  p,  which  involve  a  specific  random  variable  out  of  the  set  {{.(0)}“i  increases 
with  p.  This  fact  plays  an  important  role  in  connection  with  the  finite  dimensional 
Polynomial  Chaoses  to  be  introduced  in  the  sequel.  Furthermore,  since  random 
variables  are  themselves  functions,  it  becomes  clear  that  Polynomial  Chaoses  are 
functions  of  functions  and  are  therefore  functionals. 

The  set  of  Polynomial  Chaoses  is  a  linear  subspace  of  the  space  of  square- 
integrable  random  variables  ©,  and  is  a  ring  with  respect  to  the  functional  mul¬ 
tiplication  r,r,(w)  =  rp(w)r,(w).  In  this  context,  square  integrability  must  be 
construed  to  be  with  respect  to  the  probability  measure  defining  the  random  vari¬ 
ables.  Denoting  the  Hilbert  space  spanned  by  the  set  {£,(0)}  by  ©((),  the  resulting 
ring  is  denoted  by  $0((),  and  is  called  the  ring  of  functions  generated  by  ©(()• 
Then,  it  can  be  shown  that  under  some  general  conditions,  the  ring  $ e(()  is  dense 
in  the  space  ©  (Kakutani,  1961).  This  means  that  any  square-integrable  random 
function  (0  — *  R)  can  be  approximated  as  closely  as  desired  by  elements  from  $©(()■ 
Thus,  any  element  p(0)  from  the  space  ©  admits  the  following  representation, 

M)  -  £  £  £  «;':•£  r,(u*) . U«)),  (23) 

S>0  n,+...+n,=j  f, . f, 

where  r,(.)  is  the  Polynomial  Chaos  of  order  p.  The  superscript  n,  refers  to  the 
number  of  occurrences  of  {,.(0)  in  the  argument  list  for  T,(.).  Also,  the  double 
subscript  provides  for  the  possibility  of  repeated  arguments  in  the  argument  list  of 
the  Polynomial  Chaoses,  thus  preserving  the  generality  of  the  representation  given 
by  equation  (23).  Briefly  stated,  the  Polynomial  Chaos  appearing  in  equation  (23) 
involves  r  distinct  random  variables  out  of  the  set  {fi(0)}“,,  with  the  k's  random 
variable  {*(0)  having  multiplicity  n*,  and  such  that  the  total  number  of  random 
variables  involved  is  equal  to  the  order  p  of  the  Polynomial  Chaos.  The  Polynomial 
Chaoses  of  any  order  will  be  assumed  to  be  symmetric  with  respect  to  their  argu¬ 
ments.  Such  a  symmetrization  is  always  possible.  Indeed,  a  symmetric  polynomial 
can  be  obtained  from  a  non-symmetric  one  by  taking  the  average  of  the  polynomial 
over  all  permutations  of  its  arguments.  The  form  of  the  coefficients  appearing  in 
equation  (23)  can  then  be  simplified,  resulting  in  the  following  expanded  expression 
for  the  representation  of  random  variables, 

p(0)=  a0  T0  +  ]£  a<i^i({i)(®)) 

■  i=i 


(24) 
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+  E  E  *..,r»te,w.u«)) 

l|Sl  tj  =  J 
oo  M  *1 

+  £  £  E  (.,(«), &M) 

M  =  1  »?  =  l  «j=sl 

^  52  52  52  52 

>|Xl  i]sl  ij si  »«  si 

where  Tp(.)  are  successive  Polynomial  Chaoses  of  their  arguments,  the  expansion 
being  convergent  in  the  mean-square  sense.  The  upper  limits  on  the  summations 
in  equation  (24)  reflect  the  symmetry  of  the  Polynomial  Chaoses  with  respect  to 
their  arguments,  as  discussed  above.  The  Polynomial  Chaoses  of  order  greater  than 
one  have  mean  zero.  Polynomials  of  different  order  are  orthogonal  to  each  other;  so 
are  same  order  polynomials  with  different  argument  list.  At  times  in  the  ensuing 
developments,  it  will  prove  notationally  expedient  to  rewrite  equation  (24)  in  the 
form 

mW  -  £  a>*vU(0)]>  (25) 

j=o 

where  there  is  a  one-to-one  correspondence  between  the  functionals  <!»[.]  and  Tf.], 
and  also  between  the  coefficients  a,  and  a, appearing  in  equation  (24).  Implicit 
in  equation  (24)  is  the  assumption  that  the  expansion  (24)  is  carried  out  in  the  order 
indicated  by  that  equation.  In  other  words,  the  contribution  of  polynomials  of  lower 
order  is  accounted  for  first. 

Throughout  the  previous  theoretical  development,  the  symbol  8  has  been  used 
to  emphasize  the  random  character  of  the  quantities  involved.  It  will  be  deleted 
10  the  ensuing  development  whenever  the  random  nature  of  a  certain  quantity  is 
obvious  from  the  context. 

As  defined  above,  each  Polynomial  Chaos  is  a  function  of  the  infinite  set  {(,} , 
's  therefore  an  infinite  dimensional  polynomial.  In  a  computational  setting, 
however,  this  infinite  set  has  to  be  replaced  by  a  finite  one.  In  view  of  that,  it 
*«ems  logical  to  introduce  the  concept  of  a  finite  dimensional  Polynomial  Chaos, 
pecifically,  the  n-dimensional  Polynomial  Chaos  of  order  p  is  the  subset  of  the 
olynotnial  Chaos  of  order  p,  as  defined  above,  which  is  a  function  of  only  n  of  the 
^correlated  random  variables  As  n  — »  co,  the  Polynomial  Chaos  as  defined 
Previously  is  recovered.  Obviously,  the  convergence  properties  of  a  representation 
^ed  on  the  n-dimensional  Polynomial  Chaoses  depend  on  n  as  well  as  on  the 
choice  of  the  subset  )<Vxi  °'jt  of  the  infinite  set.  In  the  ensuing  analysis,  this 
cfioice  will  be  based  on  the  Karhunen-Loeve  expansion  of  an  appropriate  random 
Process.  Since  the  finite  dimensional  Polynomial  Chaos  is  a  subset  of  the  (infinite- 
'mensional)  Polynomial  Chaos,  the  same  symbol  will  be  used  for  both,  with  the 
‘rnension  being  specified.  Note  that  for  this  case,  the  infinite  upper  limit  on  the 
*^mmations  in  equation  (24)  is  replaced  by  a  number  equal  to  the  dimension  of  the 
o'ynomials  involved.  For  clarity,  the  two-dimensional  counterpart  of  equation  (24) 
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is  rewritten,  in  a  fully  expanded  form,  as 

fi[0)  =  cq  F o  +  “t  r,({,)  +  aj  ri({2)  (26) 

.+  ®u  +  aurj(f3,(i)  +  aMr2(6,<a) 

+  <*in  +  <*5ii  FaUi.fi.fi)  +  “m  r3(6,6,6) 

+  “111  rate. 6. 6)  .... 

In  view  of  this  last  equation,  it  becomes  clear  that,  except  for  a  different  indexing 
convention,  the  functionals  '!'[.]  and  rj.j  are  identical.  In  this  regard,  equation  (26) 
can  be  recast  in  terms  of  'I'jj.J  as  follows 

fj(0)  =  jfl'I'o  +  “3' Vt  +  “3  +  0^4  +  ii'bi 

+  +  ar' Vi  +  as#*  +  at'ffa  +  ...  ,  (27) 

from  which  the  correspondence  between  '!'{.)  and  T[.}  is  evident.  For  example,  the 
term  ajuTatetete)  of  equation  (26)  is  identified  with  the  term  a*?' fr?  of  equation 
(27). 


3.2.2  Construction  of  the  Polynomial  Chaos 

A  direct  approach  to  construct  the  successive  Polynomial  Chaoses  is  to  start  with 
the  set  of  homogeneous  polynomials  in  {{,(£)}  and  to  proceed,  through  a  sequence 
of  crthogonalization  procedures.  The  zeroth  order  polynomial  is  a  constant  and  it 
can  be  chosen  to  be  1.  That  is 

r0  =  1  .  (28) 

The  first  order  polynomial  has  to  be  chosen  so  that  it  is  orthogonal  to  all  zeroth 
order  polynomials.  In  this  context,  orthogonality  is  understood  to  be  with  respect  to 
the  inner-product  defined  by  equation  (2).  Since  the  set  {£,}  consists  of  zero-meao 
elements,  the  orthogonality  condition  implies 

rite)  =  (i  • 

The  second  order  Polynomial  Chaos  consists  of  second  order  polynomials  in  {£,}  that 
are  orthogonal  to  both  constants  and  first  order  polynomials.  Formally,  a  second 
order  polynomial  can  be  written  as 

—  “0  +  “i,  £1,  + 

where  the  constants  are  so  chosen  as  to  satisfy  the  orthogonality  conditions.  The 
second  of  these  requires  th«* 

<r»te, .<.,)&>  -  0.  (3D 

This  leads  to  the  following  equation 

+  Oijdijij  =  0  . 


(32) 
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Allowing  h  to  be  equal  to  i,  and  i,  successively,  permits  the  evaluation 
coefficients  a,,  and  a,,  as 

o„  =  0  ,  a„  =  0 

The  first  orthogonality  condition  yields 

uo  +  =  0  . 

Equation  (34)  can  be  normalized  by  requiring  that 

a,, i,  —  1  - 


This  leads  to 


flo  —  **  ^i|«»  • 

Thus,  the  second  Polynomial  Chaos  can  be  expressed  as 


of  the 

(33) 

(34) 

(35) 

(36) 


- 


(37) 


Id  a  similar  manner,  the  third  order  Polynomial  Chaos  has  the  general  form 

=°0  +  +  *h(h  +  «*.&  +  *.■»&&» 


(38) 


with  conditions  of  being  orthogonal  to  all  constants,  first  order  polynomials,  and 
second  order  polynomials.  The  first  of  these  conditions  implies  that 


=  o  • 

(39) 

That  is, 

(40) 

<*o 

The  second  condition 

implies  that 

o 

II 

A 

v> 

F* 

V/ 

e» 

V 

V 

rt 

u, 

V 

(41) 

which  leads  to 

+  0, +  a'l'Vj  <-^'l  (" >  ‘ 

(42) 

The  last  orthogonality  condition  is  equivalent  to 

<r3((,. ,  <f<3  •  £i3  )  (’!'*  —  3 , 

(43) 

*hich  gives 


°0  ^i,i|  "b 

+  6’»  ^  ~ 


(44) 


u 


The  above  equations  can  be  normalized  by  requiring  that 

Then  equation  (42)  becomes 

d-  4-  +  <£<,£, i£<j£»4 =  0 


Substituting  for  the  expectations  in  equations  (46)  and  (44)  yields 

0  , 


and 


no  6,,;,  +  a, 

4-  a 

+  a 


4- 

4-  0,16,,., 

4- 

4"  6,,,,  5,,i, 

4- 

6<i<.6,,i,  — 

<i*i  1 

6<1<16<4<1 

4-  6jl;,6j,„ 

4- 

6<1<16<J<4  ) 

«i*i  t 

6<1<16<4<1 

4-  6,4,4 6,ltl 

4- 

6<1<16<1<4  ) 

•j'>  1 

6<3<16<4<1 

4-  6,,,, 6,1,1 

4- 

6<3<16<1<4  ]  — 

=  0. 


(45) 

(46) 


Due  to  the  Gaussian  property  of  the  set  {(, },  the  following  equation  holds 

~  4"  ^tii)6i)t(  4  »« i s  *  (47) 


(48) 


(49) 


Substituting  for  a0  from  equation  (40),  equation  (49)  can  be  rewritten  as 

a'i >j  [  6i|it^,jit  4'6.,iJ^t,(<  )  4-  a,, i3  [  4-  ) 

4-a,-,,-,  [  4-  £.1,1<5.,u  )  =  0  .  (50) 

From  equation  (50),  the  coefficients  a,,,-,,  a,-,,,,  and  a,-,,-,  can  be  evaluated  as 


°>i,i  =  0 

a,,ij  —  0 

a.ji]  =  0  . 


Using  equation  (49)  again,  it  is  found  that 

o0  =  0  . 


(5i; 


(52) 


Equation  (48)  can  be  rewritten  as 

6i|„(o,,  4- 4-  6,,;, (a;,  +  6,JI})  +  <5,,„(a,i  4-  6.,,,)  =  0,  (53) 


/  j 


from  which  the  coefficients  a,,,  a„,  and  a„  are  found  to  be. 


(54) 


The  third  order  Polynomial  Chaos  can  then  be  written  as 

*(•*>£<»)  =  £>i  4ii  (•!  -  ^I'i  —  (■>  ~  £*j  ^'i*i  ■ 


(55) 


After  laborious  algebraic  manipulations,  the  fourth  order  Polynomial  Chaos  can  be 
expressed  as 


>  ('»  (<»•(•* ) 


=  (i|  (it  (it  ('• 

—  (it  (it  ^ilit  ( 1 1  (ij  ^'J't  ('I  ^*1'1 

(ij  Sit  i,  —  (it  (i,  (t,it  ~  (<J  (lt  ^'1  '3 

+  ^'1  '3 4"  (iiii^iVt  4"  * 


(56) 


It  is  readily  seen  that,  in  general,  the  n'*  order  Polynomial  Chaos  can  be  written  as 


£  M)r 

*wtn 


£  n  <  n  (i< > 

to, . ;»)  *=i  <=•■+« 


n  even 

(57) 

£  (~ i)r-'  n  &  <  n  &>> 

»(■! . i«l  *=l  txr+l 


n  odd 


where  *•(.)  denotes  a  permutation  of  its  arguments,  and  the  summation  is  over  all 
*uch  permutations  such  that  the  sets  6, }  is  modified  by  the  permutation. 

Note  that  the  Polynomial  Chaoses  as  obtained  in  equations  (28),  (29),  (37),  (55) 
(56)  are  orthogonal  with  respect  to  the  Gaussian  probability  measure,  which 
^akes  them  identical  with  the  corresponding  multidimensional  Hermite  polynomials 
(Grad,  1949).  These  polynomials  have  been  used  extensively  in  relation  to  problems 
m  turbulence  theory  (Imamura  et.al,  1965a-b).  This  equivalence  is  implied  by  the 
orthogonality  of  the  Polynomial  Chaoses  with  respect  to  the  inner  product  defined 
equation  (2)  where  dP  is  the  Gaussian  measure  d£,  where  t  denotes  the 


vector  of  n  random  variables  {(„ }  J_,  Tliis  measure  is  exactly  the  weighing  function 
with  respect  to  which  the  Hermite  polynomials  are  orthogonal  in  the  Lj  sense  (Oden, 
1979).  This  fact  suggests  another  method  for  constructing  the  Polynomial  Chaoses, 
namely  from -the  generating  function  of  the  Hermite  polynomials.  Specifically,  the 
Polynomial  Chaos  of  order  n  can  be  obtained  as 

r.Ki . =  Mr  (58) 

The  first  two  terms  in  equation  (25)  represent  the  Gaussian  component  of  the 
function  fi{0).  Therefore,  for  a  Gaussian  process,  this  expansion  reduces  to  a  single 
summation,  the  coefficients  a,,  being  the  coefficients  in  the  Karhunen-Loeve  expan¬ 
sion  of  the  process.  Note  that  equation  (25)  is  a  convergent  series  representation  for 
the  functional  operator  A(.]  appearing  in  equation  (23).  For  a  given  non-Gaussian 
process  defined  by  its  probability  distribution  function,  a  representation  in  the  form 
given  by  equation  (25)  can  be  obtained  by  projecting  the  process  on  the  successive 
Homogeneous  Chaoses.  This  can  be  achieved  by  using  the  inner  product  defined 
by  equation  (2)  to  determine  the  requisite  coefficients.  This  concept  has  been 
successfully  applied  in  devising  efficient  variance  reduction  techniques  to  be  coupled 
with  the  Monte  Carlo  simulation  method  (Chorin,  1971;  Maltz  and  Hitzl,  1979). 


4.  PROJECTION  ON  THE  HOMOGENEOUS  CHAOS 

In  this  section  the  Karhunen-Loeve  expansion  and  the  Polynomial  Chaos  expansion 
presented  earlier  are  implemented  into  a  stochastic  finite  element  method  which 
features  a  number  of  similarities  with  the  deterministic  finite  element  method. 
Specifically,  the  geometric  interpretation  of  the  finite  element  method  as  a  projection 
in  function  space  is  preserved. 

Equation  (6)  constitute  the  starting  point.  Assuming  that 

n(x,w)[.j  =  o(x,0)  R(x)[.)  ,  (59) 

and  expanding  ct(x,0)  in  a  Karhunen-Loeve  scries  gives 

(  L(x)  +  £  c  a„(x)  R(x)  )  (u(X,0)]  =  f(x,8)  .  (60) 

Assuming,  without  loss  of  generality,  that  m(x,0)  is  a  second  order  process,  it  lends 
itself  to  a  Karhunen-Loeve  expansion  of  the  form 

L 

U(M)  =  £  ei  XA°)  6;(x)  .  (61) 

;=i 

JD  C„  (xi  .  x2)  6j(x2)  dxj  =  t,  b,(x i)  , 


where 


(62) 
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Vj(0)  -  ~  [  u{x,0)  b,{x)  dx  . 

e,  J D 


Obviously,  the  covariance  function  Cuu{X|,x2)  of  the  response  process  is  not  known 
at  this  stage.  Thus,  c,  and  ftj(x)  arc  also  not  known.  Further,  u(x,0),  not  being 
a  Gaussian  process,  the  set  \j(0)  is  not  a  Gaussian  vector.  Therefore,  equation 
(61)  is  of  little  use  in  its  present  form.  Relying  on  the  discussion  concerning  the 
Homogeneous  Chaos,  the  second  order  random  variables  \j(0)  can  be  represented 
by  the  mean-square  convergent  expansion 

Xl(0)  =  a[?  ru  +  £;  a!;'  r.**,,) 

■i=i 

+  £  £  ,«...<„)+£  E  E«%L  raUi„Mt) 

•j“>  I J  —  I  »1  =  I  «7  =  l  Ijsl 

+  £  it  i  £  r<((:  •  ^  .  (j  +  ....  (6-i) 

»J  =  j  »>=1  ‘3-1  »«  =  1 

where  a]'*  ,,  arc  deterministic  constants  independent  of  0  and  Tp(  ( ,,  ,  ...  ,  ) 

is  the  p'*  order  Homogeneous  Chaos.  Equation  (64)  is  truncated  after  the  P'1' 
polynomial  and  is  rewritten  for  convenience,  as  discussed  in  equation  (27),  in  the 
following  form, 

X,(0)  =  £  *!'*  .  (65) 

.  =  0 

where  x[^  and  'I'.KCll  arc  identical  to  a-J1,,  and  Tp(( . („)■  respectively.  In 

equation  (65),  P  denotes  the  total  number  of  Polynomial  Chaoses  used  in  the 
expansion,  excluding  the  zeroth  order  term.  Given  the  number  Af  of  terms  used 
in  the  Karh  uncn-Locvc  expansion,  and  the  order  p  of  Homogeneous  Chaos  used,  P 
may  be  determined  by  the  equation 

p  =  1  +  £-t  ft  •  (6<>) 

J=1  4'  r=0 

Substituting  equation  (G5)  for  XjWi  equation  (61)  becomes 

«(M)  =  EE  x';)  *.[{(rllcJ-(x),  (67) 

;=!  t—0 


C;(X)  =  tj  bj(x)  . 

Changing  the  order  of  summation  in  equation  (67)  gives 

u(x,0)  =  £  *,|Ur}}  E  cj(*) 

«2=0  J=  I 


!<• 


=  £  +.lttUK(*) 


d,(x)  =  £  o(x)  . 

k=t 

Substituting  equation  (69)  for  u(x,0),  equation  (GO)  becomes 


(  L(x)  +  £  <„  «in(x)  R(x)  )  £ 

V  n=!  /  ;=0 


X)  =  /(X),  (71) 


where  reference  to  the  parameter  3  was  eliminated  for  notalional  simplicity. 

The  response  u(x,0)  can  be  completely  determined  once  the  functions  d,(x)  are 
known.  In  terms  of  the  eigenfunctions  6;(x)  of  the  covariance  function  of  u(x,0), 
d,(x)  can  be  expressed  as 

d,(x)  =  52  x>’]  eJ  Mx) 

j=i 

=  £  ylJ)  Mx)  .  (72) 

J=> 

Equation  (71)  may  be  written  in  an  alternative  form 

52  *-,[{*,})  L(X)  (,/;(x)j  +  52  £  6  *iltfr)i  R(x)  [d,(x)|  =  /(x).  (73) 

}~0 

This  form  of  the  equation  shows  that  d,(x)  belongs  to  the  intersection  of  the  domains 
of  R(x)[.j  and  L(x)(.|.  Then,  following  the  standard  deterministic  finite  element 
method,  the  function  d,(x)  may  be  expanded  in  an  appropriate  function  space  as 

<*;(x)  =  £  dk,  gk{x)  .  (74) 

fc=l 

Then,  equation  (73)  becomes 

£  £  *>,  L(x)  |*(x)]  (75) 

J=0  k=\ 

+  £  £  m  £  dkl  R(x)  Isr*(x)J  =  /(x)  . 

;=0  1=1  k=l 

Equation  (75)  may  be  rearranged  to  give 

£  £  dk)  (  *,({{,}]  L(x)  lst(x)] 

;=0  1=1 


+  £  fcW  *i(Kr)l  R(x)  ls*(x)l|  =  /(x). 

»=1 


7/ 


Multiplying  both  sides  of  equation  (76)  by  gi(x)  and  integrating  throughout  yields 

E  E  dki  f  L  L(x)  [#fc(x)|  <?i(x)  dx 

j=Q  i=t  L  J U 

+  E  £*  L  R(*)  U*(x)i  s((x) 

.=sl 

=  JD  /(*)  9i(x)  dx  ,  l  —  1  N  . 

L“  =  Id  k*(x)l  9‘^  dx 


Setting 


R'*'  =  Id  ^(x)!  ffi(x)  a,(x!  dx 


~  In  9'^  dx  ' 


equation  (77)  becomes 


p  N 

EE 

j= o  k=l 


M 


♦jHMI  u,  +  £  (.{(?)  *,•!{&}]  r.-», 


=  f<  ,  /=  l . N  . 


(77) 

(78) 


(79) 


(80) 


(81) 


Note  that  the  index  j  spans  the  number  of  Polynomial  Chaoses  used,  while  the 
index  k  spans  the  number  of  basis  vectors  used  in  Cm.  Multiplying  equation  (81) 
*>y  4'm[{{,}),  averaging  throughout  and  noting  that 

<*>1{C)1  =  Sjm  <'l'U{?r)]>  ,  (82) 

one  can  derive 

E  <*mKMl>  w*m  +  E  E dn  E  <.(0)^({Mi'M{M!>R.*< 

*=!  J=0k  =  l  i=| 

=  <fi  /  =! . N,m  =  l,...,P.  (83) 

Introducing 

C,Jm  B  <{, ;#,[«,}}  #»[{{, }J>  ,  (84) 

and  assuming,  without  loss  of  generality,  that  the  Polynomial  Chaoses  have  been 
normalized,  equation  (83)  becomes 

n  p  N  m 

E  L*,  dim  +  EE  d*>  E  R-*<  C\jrn  —  ^f|  I 

4=1  J=0  4=1  1=1 

/  =  i,...,yv  ,  m  =  1 . P  . 


(85) 


For  a  large  number  of  index  combinations  the  coefficients  r,;„,  arc  identically  zero. 
Equation  (S4)  was  implemented  using  the  symbolic  manipulation  program  MAC- 
SYMA  (19S6).  Forming  equation  (83)  for  all  P  values  of  m,  produces  a  set  of  N  xP 
algebraic  equations  of  the  form 


(  G  +  R  ]  d  =  h  ,  (86) 

where  G  and  R  are  block  matrices  of  dimension  N  x  P.  Their  mj‘1'  blocks  are 
jV-dimensional  square  matrices  given  by  the  equations 

Gmj  =  6m]  L  ,  (87) 

and 

=  £  c,}m  R,  .  (88) 

1=1 

In  equations  (S7)  and  (SS),  L  and  R,  denote  (V-dimensional  square  matrices  whose 
kl,h  element  is  given  by  equations  (78)  and  (79),  respectively.  In  equation  (86),  h 
signifies  the  N  x  Af  vector  whose  m"'  block  is  given  by  the  equation 

hm  ,  <f  *„[{*,}]>  .  (89) 

The  Af-dimensional  vectors  dm  can  be  obtained  as  the  subveclors  of  the  solution 
to  the  deterministic  algebraic  problem  given  by  equation  (SC).  Once  these  coeffi¬ 
cients  are  obtained,  back  substituting  into  equation  (69)  yields  an  expression  of  the 
response  process  in  terms  of  the  Polynomial  Chaoses  of  the  form 

*,IKr)]  •  (90) 

J=0 

Based  on  equation  (90),  realizations  of  the  random  response  vector  can  be  computed 
from  realizations  of  the  random  variables  {&■}•  Also,  statistical  moments  of  the 
random  response  vector  can  be  evaluated  using  the  inner  product  defined  in  equation 
(2). 


5.  Numerical  Examples 

The  preceding  development  of  the  stochastic  finite  element  method  was  applied 
to  a  number  of  problems  from  engineering  mechanics.  The  first  step  in  the  solution 
of  any  of  these  problems  was  the  solution  of  the  eigenvalue  problem  associated  with 
the  Karhuncn-Loeve  expansion.  Following  that,  the  coefficients  in  the  Polynomial 
Chaos  expansion  for  the  solution  process  were  computed.  Finally,  various  statistics, 
as  well  as  the  probability  distribution  of  the  solution  process  were  numerically 
evaluated.  Figure  (1)  shows  a  thin  plate  whose  modulus  of  elasticity  is  assumed  to  be 
a  two-dimensional  random  process.  The  plate  is  analyzed  using  the  stochastic  finite 
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Figure  1:  Plate  with  Random  Rigidity;  Exponential  Covariance  Model. 

element  formulation  described  above.  Figure  (2)  compares  some  of  the  coefficients 
inequation  (90)  for  various  levels  of  approximation;  note  the  excellent  convergence. 
Figure  (3)  shows  the  variation  of  the  standard  deviation  of  the  response  against  the 
standard  deviation  of  the  material  property  again  for  various  levels  of  approximation. 
Finally,  figure  (4)  shows  the  probability  distribution  of  the  response  variable  at  the 
free  corner  of  the  plate. 


6.  Conclusions 

A  method  for  the  solution  of  differential  equations  with  random  orocesses  as  co¬ 
efficients  was  discussed.  The  method  relics  on  viewing  the  camio  n  aspect  of  the 
problem  as  an  added  dimension,  and  on  treating  random  variables  and  r.TO^esses 
«  functions  defined  over  that  dimension.  In  this  manner,  a  formulation  for  the 
stochastic  finite  clement  method  was  derived  which  could  be  construed  as  a  natural 
extension  of  the  deterministic  finite  element  method.  Finite  element  representation 
along  the  random  dimension  was  achieved  via  two  spectral  expa  isions.  One  of  them 
was  used  to  represent  the  coefficients  of  the  differential  equation  which  model  the 
random  material  properties,  the  other  was  used  to  represent  the  random  solution 
Process.  The  new  concepts  were  implemented  using  a  number  of  computational 
models  for  simple  engineering  systems.  The  convergence  of  the  discussed  approx¬ 
imations  was  demonstrated  numerically.  Probability  distribution  functions  of  the 
response  variables  were  obtained. 

The  present  formulation  can  be  viewed  as  a  definite  step  towards  a  unification 
of  various  finite  element  techniques.  Indeed  it  consists  of  generalizing  the  concepts 
of  finite  element  approximation  to  abstract  spaces,  of  which  the  usual  euclidian 
’pace  is  a  special  case.  The  deterministic  case  can  then  be  regarded  as  a  digress, on 
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Figure  2:  Linear  Interpolation  of  tlie  Nodal  Values  of  the  Vector  c,  of  Equation 
(5.73)  for  the  Rectangular  Plate  Stretching  Problem,  i  =  2;  Longitudinal  Displacement 
Representation;  2  Terms  in  K-L  Expansion,  M  =  2;  P  =  3,0,  JO. 
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Figure  3:  Normalized  Standard  Deviation  of  Longitudinal  Displacement  at  Corner  A  of 
the  Rectangular  Plate,  versus  Standard  Deviation  of  the  Modulus  of  Elasticity;  amtz  s 
0.433;  Exponential  Covariance;  Polynomial  Chaos  Solution 
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DISPLACEMENT 

Figure  4:  Longitudinal  Displacement  at  the  Free  End  of  tlic  Rectangular  Plate; 
Probability  Density  Function  Using  30.000-Sainple  MSC,  and  Using  Third  Order 
Homogeneous  Cliaos;  Four  Terms  in  the  K-L  Expansion;  Exponential  Covariance. 

of  this  formalism  to  the  particular  instance  when  the  space  of  elementary  events 
consists  of  a  single  element,  and  where  the  probability  density  function  induced  on 
the  associated  cr-algebra  is  the  uniform  distribution. 
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Pseudo  -  Simulation  Method  for  Stochastic  Problems 

B.A.  Zeldin1 2,  P.D.  Spunos1 


Abstract 

A  new  numerical  method  for  problems  of  stochastic  mechanics  and  other  areas  involving  a  small 
number  of  random  parameter!  is  presented.  It  is  analogous  to  the  Monte-Carlo  simulation  method  and 
quite  more  efficient.  As  an  example,  the  eigenvalue  problem  of  a  clamped -damped  beam  with  random 
rigidity  is  considered. 

Key  words:  Mome-Cario  method.  Finite  Element  method.  Finite  Difference  method,  partition, 
sample  space,  eigenvalue  problem,  mean  value,  standard  deviation. 

Introduction 

The  Monte-Carlo  simulation  method  has  been  widdy  used  in  the  area  of  stochastic  mechanics  and 
others  fields,  primarily  because  of  its  versatility  This  method  treats  a  stochastic  problem  as  an  ensemble 
of  deterministic  ones.  After  solving  a  number  of  detemumsuc  problems,  statistical  analyst!  11  performed 
10  deduce  the  necessary  parameters  charactenzing  the  random  nature  of  the  solution.  Often  this  method 
u  d»  only  option  available  to  tolve  complicated  stochastic  problems.  However,  indiscriminate  use  of 
die  method  can  not  be  advocated  due  to  its  considerable  computational  cost.  In  this  paper  a  new  numer¬ 
ical  method  for  stochastic  problems  is  presented.  In  essence,  it  is  a  Mome-Cario  simulation  method  uu- 
luing  a  limited  number  of  random  variables. 

Formulation 

Consider  a  problem  governed  by  the  equation 


L($)  u -/(?>.  <I> 

where  5  -  ...?w)  is  a  random  vector.  L(5>  is  a  mathematical  operator  desenbing  the  perfor¬ 

mance  of  this  system  which  depends  on  4-  Further.  /(§)  depends  on  the  same  parameter  set  and 
describes  the  load.  The  number  M  is  assumed  to  be  small,  and  are  statistically  independent 

random  variables. 

Solving  equation  (1)  is  equivalent  to  finding  some  Funcuon  u  «•  u(p  which  satisfies  this  equa¬ 
tion.  That  is.  for  every  realization  4  «  (?|.4j.  ••  •{«<)  of  the  random  vector  §  there  exisla  a  determinis¬ 
tic  funcuon  u  (5)  which  satisfies  the  equation  (l ). 

Consider  the  space  ft"  of  4,  x  43  x  . .  ,$M  as  the  sample  space  Cl  This  space  can  be  divided  into  N 
mbdomains  fl;,  i=l  ,...N  having  the  shape  of  M-dimensional  disjoint  rectangles  with  presented  proba¬ 
bility  mass  Next  introduce  the  set  of  funcuons  9,,  i=l,...N  such  that 


,1.  if  Sen, 

0,  otherwise 


dearly,  since  (4)  and  q^.(4)  have  disjoint  support,  unless  » *j. 


(2) 


!?(S)9,<S)9,<S>/’,(4>‘*?  “ 0  ,f  '*!•  <3) 

where  p{  ($)  is  the  probability  density  funcuon  of  4 .  and  q  (4)  is  an  arbitrary  random  v unable.  This 
meins  thar 'the  set  {9,t  J'.  ,  it  orthogonal  with  arbitrary  weight  basis  for  the  class  of  random  vanables 
wtuch  are  constant  for  every  Q,. 
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PROBABILISTIC  MECHANICS 


The  solution  of  equation  (1)  can  be  approximated  by  the  senes 

N 

u(5)  -  £c, «>,(!;)  •  (*) 

■  - 1 

where  u($)  is  assumed  to  be  an  adequately  smooth  function  of  Next,  the  scalar  product  of  two  ian- 
dom  vansbtes  a  and  b  can  be  introduced  by  the  equation 

E  [ah]  «*  (a,  b).  (3) 

Then,  the  solution  given  by  equation  (4)  can  be  constructed  as  a  projection  of  the  exact  solution  into  the 
space  Sp  { ipj}  *  .  Thus,  the  foliowing  sequence  of  equations  can  be  wntten 

p 

<2. (*)£ <7*. ($),*.($»  -  .  (6] 
i  -  1 

which,  because  of  equations  (3)  and  (5),  leads  w 

ElLWc^q))  (7) 

This  sequence  of  deterministic  equations  can  be  solved  to  determine  the  solutions  cr  Upon  determining 
Cj,  the  statistical  properties  of  the  solution  can  be  estimated  by  relying  on  equation  (4)  Specifically. 

H  ft 

£(«($)]-  £<:,/>,  and  £[«J($)1  -  £r?p,  .  (8) 

im  I  <-  1 

where  p,  -  £  Iv,  (?)  ]  -  £[*?(?)] 

tfE($)  ■  ^  .where  «  some  known  function  of  $  Mid  L<  is  linear,  equation 

(6)  leads  to  the  simple,  expression 


3pV.fr  (?)*,<?)  ] 


-  £[/(?)*,(?)]■ 


(9) 


Finite  Element  Method  *  Pseudo  Monte-Carlo  Perspective 

Examining  the  proposed  method  one  may  view  it  as  a  Finite  Element  method  for  random  media. 
That  is.  it  involves  approximations  of  the  random  variables  in  the  finite  dimensional  subspaces  defined 
by  some  finite  partition  of  the  sample  space. 

From  another  perspective,  examining  equauon  (7)  one  can  deduce  that  dvia  equation  is  equivalent 
to  the  following 


|  lM§)c(-/($))^(5)d$  -0.  (10) 

If  the  Lebesgue  integral  involved  inequation  (10)  can  be  i  nierpreted  in  the  Riemann  sense  and  all  quan¬ 
tities  in  the  above  expression  are  adequately  smooth,  the  mean-value  theorem  states  that  there  exists 
some  (,  €  0;  such  that 

Ubc,  -/(§)  HD 

This  equation  shows  that  c,  represent  just  a  solution  of  equauon  (1)  for  the  realization  |  of  the  random 
vector  §.  Then,  the  sequence  of  the  equations  (7)  can  be  interpreted  as  a  sequence  of  Mcnte-Ca/lo  sim¬ 
ulations.  In  fact,  this  sequence  of  "pseudo-simulaiions’'  is  optimal  in  the  sense  that  every  element  of  this 
sequence  represents  a  certain  region  of  the  sample  space  and  can  be  interpreted  as  the  only  outcome 
with  given  probability. 

For  the  implementation  of  this  method,  first  the  entire  sample  space  is  divided  into  some  assembly 
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of  subsets,  the  mesh  of  events,  with  given  probability  mass  Then,  certain  element*  are  chosen  from 
every  subset,  and  the  sequence  of  equauons  (I)  is  treated  as  a  detenmmsuc  one.  It  is  further  assumed 
that  these  elements  of  the  sample  space  represent  a  realization  of  the  random  vector  £ .  Then,  statistical 
analysis  is  conducted  assuming  that  the  derived  deterministic  solutions  are  the  only  possible  outcomes 
with  given  probability.  U  is  clear  that  this  aigonthm  relates  to  the  Monte- Carlo  method.  Thus,  it  should 
work  at  least  as  well  as  the  Monte-Carlo  method  Further,  it  exhibits  the  appealing  feature  of  choosing 
the  points  from  the  subsets  tn  an  optimal  way,  namely  in  the  sense  of  the  projection  represented  by  the 
equauon  (6). 

Example:  Eigenvalue  Problem  for  a  Beam 

The  proposed  method  was  applied  to  the  eigenvalue  problem  of  a  clamped-clamped  beam  of  unit 
length;  its  ngtdity  is  a  truncated  normal  random  process  with  mean  equal  to  1  and  autocorrelation  func- 
uon 


(x,-x.)J 

ffs/Cr,. x3)  -  a1  exp ( - — ■  -) . 

The  case  with  a  -  0.3  and  c=Q5  was  taken.  The  corresponding  equation  is 

(£/(x)iO*  -  Xu.  (12) 

with  boundary  condition*  u(0)  -  u(l)  -  u'(0)  -  u'(l)  -  0.. 

This  land  of  problem  is  quite  difficult  either  for  an  analytical  or  for  a  numencal  treatment  Only  a 
few  papers  are  available  on  this  topic.  A  desenpuon  of  pertinent  analytical  methods  was  presented  by 
Bruce(l968).  In  the  papers  of  Hasselman  and  Han  (1972),  and  Gngonu  (1991)  some  numencal  exam¬ 
ples  of  solution  of  stochastic  eigenvalue  problems  can  be  found  However,  these  algorithms  can  be 
ippiied  in  the  case  of  small  randomness  only  and  can  be  computationally  costly. 

In  implementing  the  proposed  method,  first  the  discretization  procedure  is  applied  to  obtain  a  finite 
dimensional  problem.  For  this  purpose  the  Finite  Difference  method  see  also  Spanos  and  Zeldin  (1992). 
is  used  This  leads  to  the  equauon 

An  -  ku.  (13) 

where  A  is  a  matrix  with  random  vanables  as  element!  and  u  is  a  random  vector.  Subsequently  the  Kat- 
hunen-Loeve  expansion,  see  also  Ghanem  and  Spanos  (1988.1991).  can  be  applied  to  represent  the 
mitnx  A  in  die  form 


A  -  A,,*?,*,  +•  +  ....  (14) 

where  ...  are  suustically  independent  random  vanables.  The  sene*  in  equauon  (14)  can 

be  truncated  beyond  order  M,  and  the  random  vector  f;  -  ({,.  ...^J  can  be  introduced  Next,  the 

random  domain  is  divided  into  a  set  of  rectangle* 

5e  i*l...  M.  j=I  ...N) 

of  presenbed  probability  mas*.  Then,  the  basis  ( p, )  f  can  be  constmcted  to  conform  with  equauon 
(2)  Next,  u  and  k  can  be  taken  in  the  form 


“(?)  "  “  £r,P,(?) 


05) 


Substituung  equation*  (14)  and  (15)  into  equation  (13).  multiplying  it  by  p  (?) ,  and  taking  the  nuuhe- 
maucal  expectation  of  the  result  yields 


where 


(A0>?‘A,  +  4^A„)  v,  a  c,vt  ,  i=l... .AT, 

??  -  andp,  -  £*,«)!>,($)<<$  - 


(16) 


Finally,  stansucal  analysis  can  be  performed  in  conjuncuon  with  equation  (8)  to  estimate  (he  moments 
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of  the  first  two  eigenvalues  and  eigenvectors 
Numerical  Results  •  Concluding  Remarks 

The  theoretical  values  of  the  first  two  eigenvalues  for  the  deterministic  case,  when  the  beam  ngidity 
is  set  equal  to  the  mean  rigidity  of  the  problem  under  consideration,  are  Xf"  ■  22.37  and 
X^"  ■  61.67.  The  corresponding  deterministic  finite  difference  approximation  with  41  node  points 
gives  X^**/  -  22,32  and  Xf*A/  -  61.33.  It  was  found  that  despite  the  considerable  variability  in  the 
rigidity  of  the  beam,  the  variability  in  the  eigenvectors  is  negligible.  However,  the  variability  in  the 
eigenvalues  is  essential  and  for  a  •  0.3 .  as  in  this  example,  is  of  the  order  of  10%.  see  Figure  1  The 
influence  of  different  was  studied.  It  was  found  that  for  this  problem  the  third  term  in  equation  (16)  is 
almost  negligible.  The  results  in  terms  of  convergence  of  this  method  for  different  order  of  partition  of 
axis  and  that  is  for  different  numbers  of  pseudo-simulations,  are  plotted  in  Figure  1(a)  for  the 
standard  deviation  of  X, ,  and  in  Figure  1(b)  for  the  standard  deviation  of  Xj.  It  is  teen  that  this  method 
gives  quite  good  approximations  even  when  the  number  of  pseudo-simulations  ia  very  small,  whereas 
the  Monte-Carlo  method  yields  reliable  results  only  if  the  number  of  the  simulations  used  is  large.  The 
mean  of  the  lint  and  second  eigenvalues  was  found  to  be  21.98  and  60.40.  respectively.  Thus,  the  mean 
values  (lightly  decline  to  a  mailer  values  compared  to  the  deterministic  case. 

It  appears  that  the  discussed  method  can  be  applied  far  treating  a  wide  class  of  problem  dealing 
with  random  variables  and  itochasac  processes.  Thus,  further  research  might  be  warranted. 
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Stochastic  Mixed  Finite  Difference  Method 

P.D-Spanos1 2,  B_A. Zeldin1 


Abstract 

Some  aspects  of  numerical  solutions  of  stochastic  mechanics  problems  are  considered.  The  Finite  Differ¬ 
ence  method  for  discretization  of  stochastic  continuous  media  problems  it  discussed  The  Neumann  expan¬ 
sion  and  perturbation  methods  used  for  solving  the  associated  system  of  algebraic  equations  are  analyzed. 
Their  dependence  on  the  mesh  size  of  the  discretisation  is  investigated,  ft  is  shown  that  a  mixed  formulation 
using  both  strain  and  stress  as  independent  vsrisbles  improves  the  performance  of  these  methods  and  reduces 
their  dependence  on  the  mesh  size. 

Key  words:  random  variable,  random  process.  Finite  Element  method.  Finite  Difference  method,  beam 
equation,  perrubation  method,  Neumann  expansion,  convergence. 

Introduction 

Recently,  considerable  attention  has  been  given  to  the  solution  of  engineering  problems  with  randomness 
in  the  spatial  domain.  Clearly,  the  special  dependence  hinders  solving  these  problems  exactly.  Thus,  approxi¬ 
mate  methods  of  solution  have  been  pursued.  The  main  idea  has  been  to  implement  the  well  known  and 
widely  used  Finite  Element  (FEM)  or  Finite  Difference  (FDM)  methods.  A  number  of  papers  have  appeared 
recently  on  (his  topic  addressing  linear  or  nonlinear  and  static  or  dynamic  problems,  see  for  example  Ghanem 
and  Spanos(1988,1991),  Liu  W.K.,etal.(1986),  Vanmarite  and  Grigonu(l983).  Takada  (1990),  Yamazaki, 
et.al-(I986). 

Usually  such  an  analysis  involves  two  steps.  Hie  first  step  is  to  introduce  some  discretization  of  the  con¬ 
tinuous  medium  problem,  which  leads  to  a  system  of  stochssuc  algebraic  equations;  the  second  step  addresses 
the  solution  of  this  system. 

In  this  paper  some  new  aspects  of  solving  continuous  media  problems  with  randomness  in  the  spatial 
domain  using  FDM  are  considered. 

Discretization  Technique* 

Several  different  techniques  for  random  problem  discretization  have  been  developed.  Typically  they  use 
the  following  random  field  represen  taboo 

f<*)  ”  £civi(x>’  <» 

■ 

where  f  (x)  is  a  random  field,  e.  are  random  variables,  and  v;  (x)  are  basts  functions  built  upon  some  appro¬ 
priate  partition  of  the  problem  uhanem  and  Sptnos<1988,1991)  used  this  form  for  stochssuc  field  representa¬ 
tion  with  basts  functions  from  Kurhunen-Loeve  expansion.  After  such  a  discretization  in  the  random  domain, 
the  Finite  Element  procedure  was  applied.  Alternatively  Takada  (1990)  and  Deodaus,  etc.(1991)  used  the  so 
called  Weighted  Integral  method  which  is  straightforward  application  of  FEM  to  the  random  problems. 

Interestingly,  stochastic  FDM  has  not  received  much  attention,  in  the  deterministic  case  often  FDM  has 
some  advantage  over  FEM  due  lo  simplicity  in  foratulihoa  and  analysis.  The  same  can  be  said  concerning 
stochastic  problems.  A  similar  formulation  was  introduced  by  Vanmarke  and  Grigoriu  (1983)  for  the  simple 
case  of  a  statically  determinant  shear  beam. 

A  Finite  Difference  approximation  for  a  differential  operator  is  next  discussed.  Let  u  (x)  be  a  stochastic 
process.  Then,  the  mean-square  derivatives  of  u  (x)  in  some  point  x.  can  be  approximated  by  using  the  cen¬ 
tral  difference  formulae 
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where  u(  ill  random  variable  equal  to  u(x()  and  h  ii  the  mesh  si2e  of  the  partition  involving  N+l  nod es. 

Using  theae  equations  the  beam  operator  Lu(x)  «  (El  (a)  u"  (a) )  can  be  approximated  by  theequa- 
tion 

Lu  (a,)  -  -L  (E  Ij.  ,uu  2  ♦  (-2  (Elj . ,  ♦  EI^J )  Uj _ ,  ♦  (EIU ,  +  4EI,  ♦  EI(. ,)  tl, 
h 

_2(EIi  +  E,i»l)ui.l*'EI,»tu..2)  .  (3) 

where  a|  are  nodal  points,  i=0.1,...N. 

Appending  to  equation  (3)  the  corresponding  approximations  of  propel  boundary  condioona.  a  system  of 
linear  stochastic  algebraic  equations  can  be  formed.  Specifically, 

Au  -  f  ,  (4) 

where  A  is  a  mama  with  random  vanables  as  elements,  u  is  a  random  vector  of  the  solution,  and  f  is  a  vector 
of  load  which  in  general  is  random 

Solution  of  the  System  of  Algebraic  Equation 

The  solution  of  the  problem  described  by  equation  (4)  is  cnicial  us  computauonal  stochastic  mechanics 
Several  methods  have  been  developed  foe  this  purpose.  One  of  them  is  the  Monte-Carlo  simulation  method 
which  has  been  widely  used.  This  method  was  incorporated  in  FEM  by  Shinoaifca  and  his  associates  (1972), 
who  developed  an  algorithm  for  the  stochastic  field  simulation  A  recent  implementaoon  of  tlus  method 
includes  the  Neumann  expansion  of  matrices  introduced  by  Yamaaaki.et  ml  (1986)  However,  this  method 
may  demand  large  computauonal  resources.  Another  quite  common  method  relies  on  a  perturbation  expan¬ 
sion  which  provide*  sufficiently  good  results  when  the  randomness  is  not  very  large  (Tikada(I991),  Liu. 
et,al.(1990»  Otherwise  this  method  can  give  erroneous  results 

Alternatively,  the  Neumann  expansion  method  can  be  utiliied.  In  fact  that  is  some  son  of  generalization 
of  the  pemirbauon  method.  It  was  used  by  several  mvesugations  such  as  in  Ghanem  and  Spanos  (1988. 
1991),  Adomtan  and  Malkian  (1979)  and  YamuaJu.  et  al  (1988)  Let  the  matrix  A  in  the  equation  (4)  be 
expressed  as  A  «  Aq  +  A  with  Aq  -  (A)  and  A  »  A-Aq  Substituting  this  expression  into  equa- 
uon  (4)  and  tnuluplying  it  by  Aq  qives 

(I  +  A"’a)u  -  A’'f  (5) 

Further,  it  is  possible  to  introduce  the  following  equality 


“  “  L  (-‘)kfA0  A)  Ao’f  “  A0f"  A0AAo  f+  <A0  A>  Ao  f- 


Eqyguqti  (6)  has  only  a  formal  meaning  as  it  is  assumed  that  this  senes  converge  io 
(I  +  Aq  A)  AqI  .  Examinauon  of  Uus  reveals  that  the  convergence  problem  is  nor  a  trivial  one  One  of 

the  sufficient  conditions  for  this  is  the  following  almost  surely  inequality 

Ik'AlMO  <?> 

Note  that  the  condition  given  by  equsuon(7)  is  too  restrictive  and  difficult  for  analysis  Adomian  and  Malk- 
ian  (1979)  considered  linear  differenual  operator  L  “  ^»|(^)  wu8  randomness  in  ihe  lastcoelficient  Sq 


Convergence  of  iheje  senes  was  proved  and  the  enot  of  truncauon  waa  estimated  in  that  case  However,  in  a 
number  of  paper*  this  method  was  applied  even  when  other  coefficients  a  ,  ml...  N  were  random;  see  Gha¬ 
nem  and  Spanos  (1988,1991),  Adomian  (.1979),  Yamizaiu  and  ere  (1988)  It  is  not  difficult  to  note  that  for 
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continuous  problem  operates  d/dt  is  not  bounded .and  for  discrete  care  tt  behave*  a*  !/b.  where  h  u  a  mesh 
*iie.  Thus,  for  the  beam  problem  the  norm  !j  Aj  A||  may  grow  ai  1/h4. 

One  can  show  that  the  perturbatton  method  represents  the  Neumann  expansion  truncated  to  the  first  2 
terms.  Then,  tt  can  be  concluded  that  whenever  the  mesh  of  the  discretization  ts  sufficiently  small,  both  the 
perturbation  and  Neumann  expansion  methods  of  solution  of  a  system  of  algebraic  equations  may  behave 
poorly  due  to  the  unbaindednesa  of  the  differentisl  operaior. 


Mixed  Method 


It  is  possible  io  overcome  ihe  convergence  obstacle  using  the  concept  of  Mixed  FDM.  The  applicability 
of  this  concept  has  been  demonstrated  for  the  solution  of  deterministic  problems  with  non-smooth  coeffi¬ 
cients.  In  fact  this  method  allows  the  elimination  of  non-smooth  coefficients  in  the  left  hand  side  of  the  differ¬ 
ential  equations  describing  the  problem.  This  method  can  be  helpful  and  for  the  stochasuc  case  where  it  is 
necessary  to  find  an  algorithm  allowing  the  use  of  the  Neumann  expansion  procedure.  Indeed.  Tor  the  beam 
problem  the  new  variables  v,  (x)  -  u(x)  and  Vj(x)  -  EIu”  (x)  can  be  introduced.  Then  the  beam 
equauon  can  be  written  in  the  form 


v"  —  Bv  *  q. 


(8) 


where  v 


B 


and  q  ■ 


is  associated  with  a  term 
on  the  mesh  size  can  be 


El(x) I 

Lo  0 

As  one  can  see  from  equation  (10),  the  random  coefficient  of  this  i 
which  does  not  involve  differentiauon.  Therefore,  the  dependence  of  || 
avoided 

To  show  the  applicability  of  the  method  discussed,  consider  a  clamped -clamped  beam  subjected  to  deter¬ 
ministic  force  q=l  with  the  bending  rigidity  being  a  normal  random  process  as  shown  tn  Figure  1  The  bend¬ 
ing  of  this  beam  can  be  described  by  equaucxi  (S)  with  boundary  conditions  Vj  (0)  ■V|(|)i 
=v  ^  (0)  -  v'  J  (1)  m  o  To  solve  this  problem,  FDM  can  be  applied;  it  results  in  an  equation  similar  to 
equation  (4).  Consider  a  case  where 


M«)  ■  E(|x)  •  <k(x))»l  and  (It (*t)  k (x^) )  -  a2  —  -) 

where  a  »  0.3  and  c  -  1.0.  Meshes  with  10  and  20  nodes  were  chosen.  The  second  stausucal  moments  of 
displacement  and  bending  moment  were  calculated  based  on  the  expressions 

M2 

Cov(v)  -  <vv'>  -  £  <(Ag1A)*AQ/ff<Ag/((AQ1A)')  >- <u>(u>',  (9) 

iW.odd>0 


where  <u>-  £  (-Dk<(A51A)'CAglf>. 

k-0 

The  performance  of  there  stausucal  momenta  with  respect  to  the  order  of  the  truncation  was  investigated.  The 
results  obtained  foe  variance  are  plotled  in  the  Figure  2  It  is  seen  that  the  mixed  formulation  improves  the 
convergence  of  Neumann  expansion  and  removes  its  dependence  on  the  mesh  toe.  There  results  are  not 
affected  by  the  perubon  and  even  the  first  order  perturbation  method  is  acceptable.  This  is  especially  true  with 
regards  to  the  displacement  variability  even  if  the  coefficient  of  vanauon  of  the  beam  ngtdity  is  large.  How¬ 
ever,  the  bending  moment  variance  performs  poorly  compared  to  the  displacement  one  and  the  second  term  of 
the  expansion  (9)  should  be  used  It  is  seen  that  30%.  approximately,  of  the  vanebility  of  the  beam  ngtdity 
induces  10%,  approximately,  for  the  coefficient  of  variation  response.  That  coincides  with  other  results  (Gha- 
nem  and  Spanos  (1988,1991)). 
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Pigure  2.  Variance  of  the  displacement  {»  &  b)  and  bending  motneni  (c  4  d)  for  the  mesh 
10  (a  4  c)  and  20  (b  4  d)  nodes 

(A) -  Mj-2  (B)  - M2«4  (C)  -------  Mj  *  6 
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INDIRECT  SAMPLING  METHOD  FOR  STOCHASTIC  MECHANICS 

PROBLEMS 

P-D-Spanos1  B.A.Zeldin2 

INTRODUCTION 

The  Monte-Carlo  simulation  method  has  been  widely  used  in  the  field  of  stochastic  mechanics  and 
others  fields,  primarily  because  of  its  versatility.  Often  it  is  the  only  option  available  to  solve  complex 
problems.  However,  indiscriminate  use  of  the  method  can  not  be  advocated  due  to  its  considerable  com¬ 
putational  cost.  In  fact,  several  variance  reduction  techniques  have  been  developed  in  this  regard.  They 
involve  importance  sampling,  stratified  sampling,  and  others  [1], 

A  numerical  method  for  problems  of  stochastic  mechanics  and  other  areas  representing  the  solution 
by  a  small  number  of  random  parameters  is  presented.  In  essence,  it  is  a  stratified  sampling  method,  but 
more  efficient.  Alternatively,  this  new  method  can  be  viewed  as  a  Galerkin  approximation  in  the  sample 
space.  Several  examples  are  considered  involving  the  use  of  the  Loeve-Karhunen  expansion  for  stochas¬ 
tic  fields  approximation  [2,3].  The  examples  deal  with  the  evaluation  of  natural  frequencies  and  seismic 
response  of  beams  with  random  rigidity. 

FORMULATION 
Solution  Representation 

Consider  a  problem  with  random  parameters  governed  by  the  equation 

UQumftf),  (1) 


where  §  ~  (^,  f;2,  •  is  a  random  vector,  and  L(^)  is  a  mathematical  operator  describing  the 
performance  of  the  system.  Further,  /(!;)  describes  the  load,  and  £2>  •  316  statistically  inde¬ 

pendent  random  variables. 

+  Solving  equation  (1)  is  equivalent  to  finding  a  function  u  -  u  (fj)  such  that  for  every  realization 
§  —  i>  ^  2'  ■■■%  Af)  of  the  random  vector  §  there  exists  a  deterministic  function  u  )  which 
satisfies  equation  (1).  Consider  the  space  RM  of  ^  1  x  x  .  .  ^  as  the  sample  space  fl  This  space  can 
be  divided  into  N  subdomains  or  strata  {£}.,/■  1,1V }  having  the  shape  of  M-dimensional  disjoint 
rectangles  with  prescribed  probability  mass  as  shown  in  Figure  1  for  M=2.  Next,  introduce  the  set  of 
functions  or  spline  basis  [  qp(  ,  i  =  1,  ...N]  such  that 


*/($) 


1,  if  5  e  Q, 
{  ' 
0,  otherwise 


(2) 


Clearly,  <p(.(p  andcp^.(^)  have  disjoint  supports.  That  is, 

=  0  if  i  *j. 


(3) 


where  (fj)  is  the  probability  density  function  of  \ ,  q  (§)  is  an  arbitrary  random  variable.  Thus,  the 
set  { =  1,  ..  N}  is  an  orthogonal  basis  for  the  class  of  random  variables  which  are  constant  for 
every  Q..  Any  random  variable  can  be  approximated  adequately  by  the  use  of  these  basis  functions, 
provided  the  partition  of  Q  is  fine.  The  solution  of  equation  (1)  can  be  represented  as  a  linear  combina¬ 
tion  of  the  functions  {<p(.,  i  «*  1, . .  ,N\  .  That  is. 

N 

“(§)  “  Ec«*/($)-  (4) 

i « i 
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where  the  coefficients  c(-  are  to  be  determined. 

Then,  the  solution  given  by  equation  (4)  can  be  construed  as  a  projection  of  the  exact  solution  into 
the  space  spanned  by  [q>.,  i  =  1, ..  N]  .  Expressing  the  solution  in  the  form  of  the  equation  (4),  the 
induced  error  in  the  equation  (1)  can  be  made  orthogonal  to  the  space  spanned  by  { <p i  =  1 , . .  ,S]  . 
That  is,  using  the  operator  of  mathematical  expectation.  <  >.  it  can  be  found 
N 

(L  ($)  £  C.<p. ( (?)  <p . (§))>-  <f(  (?)  <py  (£) )  >  .  jm  1 . At  .  (5) 

i  -  1 


which,  because  of  equation  (3),  leads  to 

(L(§)ci<Pf(§))  =  </■(§)  <p.(§)> 


(6) 


This  sequence  of  deterministic  equations  can  be  solved  to  find  the  coefficients  c . .  Upon  deriving 
c(-,  the  statistical  properties  of  the  solution  can  be  estimated  by  relying  on  equation  (4).  Specifically, 

/v  yv 

<“($>>-  £c*pf  “d  <a2(§)>  =  .  (7) 

i  -  1  /  -  1 


where 


P,  «  <»,(§)>-  <<P?(§)  >  . 


(8) 


Similarly,  the  distribution  function  of  the  solution  can  be  found  using  the  equation 

N 

Pu(v)  =  Pr(u<v)  =  £  xv(c, )/>,-•  (9) 

/  *■  1 


where 


Xv 


(x)  =  { 


1 

0 


if  x  £  v 
if  x  >  v  ’ 


(10) 


Solution  Interpretation 

The  proposed  method  may  by  viewed  as  a  Galerkin-type  procedure  for  random  media.  Several 
authors  have  explored  the  idea  of  using  projection  procedures  in  conjuction  with  random  variables.  In 
references  [2.3.4,5,6]  this  procedure  has  been  applied  for  stochastic  mechanics  problems  with  random¬ 
ness  in  the  spatial  domain.  Due  to  the  correlation  between  the  solution  and  the  random  parameters 
describing  the  properties  of  the  structures,  this  class  of  problems  is  especially  difficult  to  solve.  In  this 
regard,  the  stochastic  field  has  been  discretized  by  the  use  of  the  Loeve-Karhunen  expansion  in  refer¬ 
ences  [2,3,4]  or  of  the  midpoint  method  in  reference  [5].  In  this  manner  the  problem  is  first  character¬ 
ized  by  a  finite  set  of  random  variables.  Then,  the  solution  can  be  derived  by  a  Galerkin  projection  into 
finite  dimensional  spaces  spanned  by  orthogonal  chaos  polynomials  as  in  references  [2,3,4],  or  just  lin¬ 
ear  functions  as  in  reference  [5].  However,  these  bases  can  yield  a  large  order  system  of  equations  which 
must  be  solved  to  determine  the  solution. 

Another  possible  basis  for  the  representation  shown  in  equation  (4)  is  given  by  equation  (2);  see 
also  reference  [6],  The  concept  of  using  spline  type  approximation  has  been  discussed  widely  in  the  area 
of  computational  mechanics  in  connection  with  the  finite  element  method.  From  this  perspective,  the 
system  of  functions  { cp . }  represents  the  simplest  spline  of  piecewise-constant  functions.  Then,  the 
proposed  method  involves  approximation  of  the  random  variables  in  the  finite  dimensional  subspace  of 
splines  defined  by  some  partition  of  the  sample  space.  Additional  advantages  of  this  representation 
relate  to  equation  (3)  since  each  term  c(-  in  the  expansion  (4)  can  be  found  independently.  Therefore, 
every  term  c(.  in  equation  (4)  can  be  readily  determined. 
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From  another  perspective,  the  use  of  piece-wise  constant  functions  makes  this  method  a  general¬ 
ized  sampling  procedure.  Indeed,  examining  equation  (6)  one  can  deduce  that  this  equation  is  equivalent 
to  the  following 

j  IL  (pc,.  -/(f)]  -  0.  (11) 

i 

If  the  Lebesgue  integral  involved  in  equation  (11)  can  be  interpreted  in  the  Riemann  sense  and  all  perti¬ 
nent  quantities  are  adequately  smooth,  the  mean-value  theorem  states  that  there  exists  some  f  €  Q . 
such  that 

LVC)ct  -ffi*).  (12) 

* 

This  equation  shows  that  c-  represents  just  a  solution  of  equation  (1)  for  the  realization  f  of  the  ran¬ 
dom  vector  f .  Then,  the  sequence  of  equations  (12)  can  be  interpreted  as  a  sequence  of  samplings.  In 
fact,  this  “indirect  sampling”  is  optimal  in  the  sense  that  every  element  of  this  sequence  represents  a  cer¬ 
tain  region  of  the  sample  space  and  can  be  interpreted  as  the  only  outcome  with  a  given  probability. 
Moreover,  as  f  £  Q.,  the  proposed  method  can  be  viewed  as  analog  to  stratified  sampling  [1],  But 
unlike  the  stratified  sampling  method  the  point  inside  every  stratum  is  computed  to  make  some  error  of 
the  approximation  of  the  given  numerical  problem  (1)  orthogonal  to  the  chosen  space  and  minimal  for  a 
given  stratification. 

To  show  this  properly  regression  analysis  can  be  applied  [7].  Any  random  variable  ©(f)  which  is 
a  function  of  f  on  the  sample  space  can  be  estimated  by  ©  using  the  set  of  random  variables 
{ <pt-,  i  =  1,  ...N)  defined  by  the  equation  (2),  where  ©  is  an  arbitrary  function  of  {tp.,  i  =  1,  ...N} 
rather  than  f .  This  estimate  provides  a  minimal  variance  for  the  difference  0  -  0 .  That  is, 

_  ~  2 

((0-0)  )  is  minimal.  (13) 

It  can  be  shown,  that  as  { tp i  =  1, . .  .N}  are  indicator  functions  of  disjoint  sets,  the  estimate  0  can  be 
mid  using  linear  regression  analysis  and  the  solution  can  be  expressed  as 

oo 

•  (14) 

i«  0 

where  0^.  =  (0|D.)  =  (©cp.)  ,  and <  I  >  denotes  conditional  expectation. 

Let  u  (f )  be  the  exact  solution  of  the  equation  (1),  and  let  «  be  an  approximation  of  this  solution. 
Define  the  error  of  such  an  approximation  by 

e-L(S)Z-U$)ifL($u-ffi).  (15) 

Then,  the  estimate  e  of  e  can  be  derived  using  equation  (14),  If  the  approximation  u  is  taken  from  a 
system  of  equations  (5),  then  e  =  0.  Thus,  the  proposed  method  ensures  that  the  error  defined  by  the 
equation  (15)  has  a  zero  mean  square  estimate  from  the  indicator  functions  of  chosen  stratification. 

Related  perspective  can  be  generated  using  some  algebra  concepts  [7].  The  vector  f  defines  a 
sigma-algebra  G  in  the  sample  space,  and  the  stratification  shown  in  Figure  1  defines  a  more  coarse 
pine  atomic  a-algebra  GjCG  with  indicator  functions  { <p(.,  /  =  1, . .  N]  .  Then,  the  above  regression 
analysis  applied  to  the  error  e  leads  to 

(e>  -  0  .  e  =  (e|G,)-0.  (16) 

In  other  words,  this  method  yields  the  minimal  error  defined  by  equation  (15)  with  respect  to  the  coarse 
o- algebra  of  given  stratification. 
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EXAMPLES 

Preliminary  Remarks 

Tlie  proposed  method  is  applied  for  die  analysis  of  die  dynamic  behavior  of  a  beam  of  unit  iengdi. 
The  beam  problem  can  be  described  by  the  equadon 

(El  (x)  u"(.\.  /) )  "  =  1/(1 •./).  (17) 


where  u  is  the  beam  deflection,  and  q  denotes  the  distributed  force  acting  on  die  beam  which  in  general 
is  taken  as  a  stochasdc  process.  The  symbol  Fl(.x)  denotes  the  beam  bending  rigidity  which  is  assumed 
to  be  a  normal  homogeneous  stochastic  process  with  mean  equal  to  1  and  autocorrelation  function 


REl(xVx2> 


(18) 


where  a  and  c  are  constants.  Thus,  randomness  is  manifested  in  this  problem  through  the  operator  and 
the  load. 

In  implementing  the  proposed  method,  first  the  approximation  of  the  stochastic  field  EI(x)  through 
a  finite  set  cf  random  variables  is  derived.  For  this  purpose,  the  Loeve-Karhunen  expansion  is  deemed 
especially  effective.  It  is  an  optimal,  in  the  mean  square  sense,  representation  of  the  field  over  the  set  of 
random  variables.  Subsequent  application  of  the  finite  difference  scheme  [8]  or  of  any  alternative  dis¬ 
cretization  scheme  leads  to  the  system  of  linear  algebraic  equations 

0VMl  +  XmAm)u  =/  .  (19) 

where  A  j,  ■  AM  are  matrices  the  dimension  of  which  depends  on  the  number  of  nodes  used  for  the 
discretization,  u  is  a  vector  representing  the  solution  at  the  nodal  points,  and  /  corresponds  to  the  force. 

Next,  specific  numerical  examples  of  application  of  the  proposed  method  are  presented;  the  numer¬ 
ical  values  ct  =  0.3  and  c  -  0.5  are  used, 
beam  Eigenvalue  Problem. 

The  eigenvalue  problem  of  a  clamped-clamped  beam  is  considered  first.  Then,  the  force  in  equation 
(17)  takes  the  form  q  (xj  =  Xu  (x) .  This  kind  of  problem  is  quite  difficult  either  for  an  analytical  or  for 
a  numerical  treatment.  Only  a  few  articles  are  available  on  tills  topic.  A  description  of  pertinent  analyti¬ 
cal  methods  was  presented  by  Boyce  [9].  In  the  papers  of  Goodwin  and  BoyceflO],  Hasselman  and  Hart 
[11]  some  numerical  examples  of  solution  of  stochastic  eigenvalue  problems  can  be  found.  Note,  that 
these  algorithms  can  only  be  applied  in  the  case  of  small  randomness  and  can  be  computationally  costly. 

In  implementing  the  proposed  method,  the  random  domain  is  divided  into  a  set  of  rectangles 
{^|  £  ij.  5  i=l  j=l,...N}  of  prescribed  equal  probability  mass.  Then,  the  basis 

{ <p(.  i  ■  1 ...  ,N]  can  be  constructed  to  conform  with  equation  (2).  Next,  u  and  X  can  be  expressed  in 
the  form 

N  N 

«($)  =  £v.<p.(S),  and  X($)  =  £c.q>.<^)  ■  (20) 

( «  1  i-l 

Substituting  equations  (20)  into  equation  (19),  multiplying  it  by  q? .  (^) ,  and  taking  the  mathemati¬ 
cal  expectation  of  the  result  yieius 

('VSl ,i4l+  mcivi  •  /V.  (21) 

where 


Finally,  statistical  analysis  can  be  performed  in  conjunction  with  equation  (7)  to  estimate  analyti 


(22) 
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cally  the  moments  of  the  first  two  eigenvalues  and  eigenvectors. 

The  theoretical  values  of  the  first  two  eigenvalues  for  the  deterministic  case,  when  the  beam  rigidity 
is  set  equal  to  the  mean  rigidity  of  the  problem  under  consideration,  are  X^1''  «*  22.37  and 
X^et  =»  61.67.  The  corresponding  deterministic  finite  difference  approximation  with  41  node  points 
gives  \f‘nd'f  =  22.32  and  =  61.33.  It  was  found  that  despite  the  considerable  variability  in  the 

rigidity  of  the  beam,  the  variability  in  the  eigenvectors  is  negligible.  However,  the  variability  in  the 
eigenvalues  is  essential  and  it  is  of  the  order  of  10%;  see  Figure  2.  The  influence  of  different  number  M 
of  used  random  variables  q .  has  been  studied.  It  was  found  that  for  this  problem  the  contribution  of  the 
terms  beyond  ^  inequation  (19)  is  negligible.  The  results  in  terms  of  convergence  of  this  method  for 
different  order  of  partition  of  axes  ^  and  that  is  for  different  values  of  number  N  or  indirect  sam¬ 
plings,  are  plotted  in  Figure  2(a)  for  the  standard  deviation  of  Xj ,  and  in  Figure  2(b)  for  the  standard 
deviation  of  X2-  It  is  seen  that  the  proposed  method  yields  quite  good  approximations  even  when  the 
value  of  number  N  is  quite  small.  However,  the  Monte-Carlo  method,  that  is  when  the  parameters  were 
sampled  arbitrarily,  yields  reliable  results  only  if  the  number  of  the  used  simulations  is  large. 

Beam  Response  to  Deterministic  Load. 

The  second  problem  involves  continuous  systems  with  random  parameters  exposed  to  deterministic 
excitation.  Specifically,  the  dynamic  response  of  a  cantilever  beam  to  earthquake-type  base  excitation  is 
considered.  In  this  case  the  force  term  in  equation  (17)  can  be  expressed  as 

q(x,t)  =  ag(t)-ii(x,  t)  -au(x,t)  ,  (23) 

where  a  is  a  coefficient  of  damping,  u  represents  the  displacement  of  the  beam  relative  to  the  base,  and 
a  (/)  is  taken  as  the  time  history  of  the  ground  acceleration  produced  by  the  North-South  component 
of  El  Centro  earthquake  recorded  on  station  No  117  and  reported  in  the  reference  [12).  It  is  shown  on 
Figure  3  for  comparisons  with  the  beam  response.  Further,  it  is  assumed  that  the  beam  has  unit  mass  per 
length.  The  discretization  of  the  beam  by  a  finite  difference  scheme  built  upon  20  nodes  in  the  spatial 
domain  in  conjuction  with  the  Loeve-Karhunen  expansion  of  the  bending  rigidity  is  used.  Then,  the 
solution  is  taken  in  the  form  of  equation  (20)  where  in  this  case  y  .  =  y(.(r)  are  deterministic  vector- 
functions.  Substituting  this  expression  into  the  resulting  equation,  multiplying  it  by  tp .  ( f) .  and  averag¬ 
ing,  an  uncoupled  system  of  deterministic  ordinary  differential  equations  is  derived.  Each  equation  of 
this  system  is  solved  numerically  using  the  central  difference  scheme.  Finally,  the  mean  value  and  the 
standard  deviation  of  the  free  end  displacement  are  determined  by  relying  on  equation  (7). 

The  time  history  of  the  free  end  displacement  of  the  cantilever  beam  having  the  mean  characteristic 
for  the  stiffness  and  damping  is  plotted  in  Figure  4(a).  Further  results  of  the  calculations  are  shown  in 
Figure  4(b,c)  for  different  value  of  number  N  of  indirect  samplings.  In  Figure  4.  Nj,  N2,  and  Na  denote 
the  number  of  strata  in  the  domain  of  ^lt  a-  respectively.  The  case  with  a  =  0.4  which  corre¬ 

sponds  to  damping  of  approximately  6%  of  critical  for  the  first  mode  of  the  system  with  deterministic 
rigidity  equal  to  the  mean  of  the  corresponding  stochastic  problem  is  considered.  Also  the  case  where  a 
is  a  random  variable  statistically  independent  from  ^  and  uniformly  distributed  between  0.1  and  0.7  is 
examined.  The  computations  show  that  3  strata  in  the  domain  of  a  can  adequately  represent  the  depen¬ 
dence  of  the  solution  on  the  damping  variability.  Also,  the  calculations  reveal  that  only  the  first  two 
components  of  the  vector  f;  influence  significantly  the  beam  response.  It  can  be  seen  that  the  dynamic 
response  of  the  beam  to  the  deterministic  excitation  is  strongly  affected  by  the  rigidity  variability. 

Beam  Response  to  Stochastic  Load. 

The  third  problem  is  described  again  by  equations  (17)  and  (23)  but  it  involves  stochastic  base  exci¬ 
tation.  Specifically,  a  (t)  is  taken  as  a  stationary  random  process.  The  proposed  indirect  sampling 
method  can  be  readily  applied  for  treating  this  problem.  The  solution  is  expressed  in  the  form  of  equa¬ 
tion  (20).  In  this  case  v  =  v.  (r)  is  not  a  deterministic  vector-function,  but  a  stochastic  vector-process. 
That  is,  v .  (f)  is  the  response  of  a  deterministic  system  to  random  excitation.  A  number  of  techniques 
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exist  for  the  solution  of  this  problem.  In  particular,  a  spectral  approach  can  be  applied  provided  that 
(t)  is  a  second  order  stochastic  process.  In  lire  latter  case  the  second-order  characteristics  of  the 
solution  can  be  determined  from  the  formulae 

N  N 

E  [ Vj  (/)  ] pi  and  ^(/j.^)  =  ]T  f?v  Op  M/’r  (24) 

i  **  1  i  ■  1 

Then,  this  approach  can  be  viewed  as  semi-analytical  method  analogous  to  the  directional  sampling.  At 
first,  the  vector  £  is  simulated,  and  then,  the  solution  for  the  given  simulation  is  calculated  using  known 
analytical  techniques  with  subsequent  application  of  some  averaging  as  in  equation  (24). 

Again,  two  cases  for  a  are  considered.  First,  a  is  a  deterministic  coefficient,  and  second  a  is  uni¬ 
formly  distributed  random  variable;  it  induces  for  each  mode  of  the  discrete  model  of  the  beam  damping 
6%,  and  2%  to  10%  of  critical,  respectively.  The  power  spectral  density  of  the  disp.  ..ement  of  the  free 
end  is  calculated  using  the  proposed  method  for  a  (/)  being  a  white  noise  process  of  unit  two  sided 
spectral  density.  The  data  for  different  numbers  or  indirect  samplings  are  plotted  in  Figure  5  together 
with  the  corresponding  solution  for  the  response  of  a  deterministic  system  with  rigidity  equal  to  the 
mean  rigidity  of  the  stochastic  system.  Figure  5  shows  that  the  randomness  of  the  system  has  a  signifi¬ 
cant  effect  on  the  system  response  variability  and  reduces  the  peaks  of  the  response  power  spectral  den¬ 
sity.  The  calculations  show  that  only  two  first  components  of  the  vector  Ij  influence  the  first  two 
moments  of  the  solution  significantly.  Further,  the  effect  of  the  damping  variability  can  be  captured 
using  only  two  strata  in  its  domain. 

Concluding  Remarks 

A  Galerkin-type  numerical  method  for  stochastic  mechanics  problems  has  been  presented.  Specifi¬ 
cally.  it  has  been  proposed  to  use  a  Galerkin  projection  into  the  space  of  simple  random  variables.  This 
space  can  be  spanned  by  the  piece-wise  constant  spline  functions  with  a  chosen  partition  of  the  sample 
space.  Further,  it  has  been  shown  that  the  proposed  method  can  be  consumed  a  generalized  sampling;  it 
is  proposed  to  call  it  indirect  sampling.  Indeed,  it  has  been  shown  that  this  method  is  closed  to  the  strat¬ 
ified  sampling  method  and  it  is  optimal  in  the  sense  of  equation  (16).  That  is.  the  approximation  of  the 
problem  from  the  space  of  simple  random  variables  produces  an  error  with  mean  and  conditional  expec¬ 
tation,  given  the  sigma-algebra  induced  by  this  partition,  equal  to  zero.  It  has  also  been  shown  tliat  this 
error  has  zero  estimate  from  die  set  of  indicator  functions  of  the  given  stratification.  Some  stcr.bastic 
mechanics  problems  have  been  studied  utilizing  die  proposed  method  in  conjuction  with  the  Loeve-Kar- 
hunen  expansion  which  is  a  versatile  tool  for  the  approximation  of  a  stochastic  field  by  a  finite  set  of 
random  variables.  These  examples  have  demonstrated  that  the  proposed  method  can  be  applied  for  treat¬ 
ing  a  broad  class  of  stochastic  mechanics  problems. 
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Figure  4.  Base  excited  vibration  of  the  top  of  the  beam:  (a)  deterministic  system 
with  mean  characteristics;  (b)  mean  value;  (c)  standard  deviation,  stochastic  system 

1) -  N  =  63(N,=9,N,=7)  2) - N  =  7  (N,=7) 

3) - N  =  45  considering  random  damping  (N]=5.  N2=3,  Na=3) 


Frequency  (sec'1) 

Figure  5.  Spectral  power  density  of  the  displacement  on  the  top  of  the  beam. 


1)  -  N  =  63  (N]=9,  N2=7)  3) - N  =  42  considering  random  damping  (N]=7.N2=3,  Nra=2) 

2)  - N  =  7  (N|=7)  4) . deterministic  system  with  mean  characteristics 
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